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Chapter 1 

Preface 

1.1 Preface to Matrix Analysis 1 



Matrix Analysis 




Figure 1.1 
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Bellman has called matrix theory 'the arithmetic of higher mathematics.' Under the influence of Bellman 
and Kalman, engineers and scientists have found in matrix theory a language for representing and analyzing 
multivariable systems. Our goal in these notes is to demonstrate the role of matrices in the modeling of 
physical systems and the power of matrix theory in the analysis and synthesis of such systems. 

Beginning with modeling of structures in static equilibrium we focus on the linear nature of the relation- 
ship between relevant state variables and express these relationships as simple matrix-vector products. For 
example, the voltage drops across the resistors in a network are linear combinations of the potentials at each 
end of each resistor. Similarly, the current through each resistor is assumed to be a linear function of the 
voltage drop across it. And, finally, at equilibrium, a linear combination (in minus out) of the currents must 
vanish at every node in the network. In short, the vector of currents is a linear transformation of the vector 
of voltage drops which is itself a linear transformation of the vector of potentials. A linear transformation 
of n numbers into m numbers is accomplished by multiplying the vector of n numbers by an m-by- n ma- 
trix. Once we have learned to spot the ubiquitous matrix- vector product we move on to the analysis of the 
resulting linear systems of equations. We accomplish this by stretching your knowledge of three-dimensional 
space. That is, we ask what does it mean that the m-by- n matrix X transforms SR™ (real n-dimensional 
space) into SR™ 1 ? We shall visualize this transformation by splitting both ffl 1 and W" each into two smaller 
spaces between which the given X behaves in very manageable ways. An understanding of this splitting of 
the ambient spaces into the so called four fundamental subspaces of X permits one to answer virtually 
every question that may arise in the study of structures in static equilibrium. 

In the second half of the notes we argue that matrix methods are equally effective in the modeling and 
analysis of dynamical systems. Although our modeling methodology adapts easily to dynamical problems 
we shall see, with respect to analysis, that rather than splitting the ambient spaces we shall be better 
served by splitting X itself. The process is analogous to decomposing a complicated signal into a sum of 
simple harmonics oscillating at the natural frequencies of the structure under investigation. For we shall see 
that (most) matrices may be written as weighted sums of matrices of very special type. The weights are 
eigenvalues, or natural frequencies, of the matrix while the component matrices are projections composed 
from simple products of eigenvectors. Our approach to the eigendecomposition of matrices requires a brief 
exposure to the beautiful field of Complex Variables. This foray has the added benefit of permitting us a 
more careful study of the Laplace Transform, another fundamental tool in the study of dynamical systems. 

-Steve Cox 



CHAPTER 1. PREFACE 



Chapter 2 

Matrix Methods for Electrical Systems 

2.1 Nerve Fibers and the Strang Quartet 1 
2.1.1 Nerve Fibers and the Strang Quartet 

We wish to confirm, by example, the prefatory claim that matrix algebra is a useful means of organizing 
(stating and solving) multivariable problems. In our first such example we investigate the response of a nerve 
fiber to a constant current stimulus. Ideally, a nerve fiber is simply a cylinder of radius a and length I that 
conducts electricity both along its length and across its lateral membrane. Though we shall, in subsequent 
chapters, delve more deeply into the biophysics, here, in our first outing, we shall stick to its purely resistive 
properties. The latter are expressed via two quantities: 

1. pi, the resistivity in ilcm of the cytoplasm that fills the cell, and 

2. p m , the resistivity in ilcm 2 of the cell's lateral membrane. 



lr This content is available online at <http://cnx.Org/content/ml0145/2.7/>. 
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A 3 compartment model of a nerve cell 




Figure 2.1 



Although current surely varies from point to point along the fiber it is hoped that these variations are 
regular enough to be captured by a multicompartment model. By that we mean that we choose a number 
N and divide the fiber into N segments each of length -J^. Denoting a segment's 

Definition 2.1: axial resistance 



Ri 



'■ N 



and 
Definition 2.2: membrane resistance 



R„ 



2ivajj 



we arrive at the lumped circuit model of Figure 2.1 (A3 compartment model of a nerve cell). For a fiber 
in culture we may assume a constant extracellular potential, e.g., zero. We accomplish this by connecting 
and grounding the extracellular nodes, see Figure 2.2 (A rudimentary circuit model). 



A rudimentary circuit model 




R 



m 



Figure 2.2 



Figure 2.2 (A rudimentary circuit model) also incorporates the exogenous disturbance, a current 
stimulus between ground and the left end of the fiber. Our immediate goal is to compute the resulting 
currents through each resistor and the potential at each of the nodes. Our long-range goal is to provide 
a modeling methodology that can be used across the engineering and science disciplines. As an aid to 
computing the desired quantities we give them names. With respect to Figure 2.3 (The fully dressed circuit 
model), we label the vector of potentials 



x\ x 2 x 3 x 4 



and the vector of currents 



V 



2/1 2/2 2/3 2/4 2/5 2/6 



We have also (arbitrarily) assigned directions to the currents as a graphical aid in the consistent application 
of the basic circuit laws. 
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R 



The fully dressed circuit model 

Ri R 




R 



m 



Figure 2.3 



We incorporate the circuit laws in a modeling methodology that takes the form of a Strang Quartet[l]: 

• (SI) Express the voltage drops via e = — (Ax). 

• (S2) Express Ohm's Law via y = Ge. 

• (S3) Express Kirchhoff's Current Law via A T y = — f. 

• (S4) Combine the above into A 1 GAx = /. 

The A in (SI) is the node-edge adjacency matrix - it encodes the network's connectivity. The G 
in (S2) is the diagonal matrix of edge conductances - it encodes the physics of the network. The / in (S3) 
is the vector of current sources - it encodes the network's stimuli. The culminating A T GA in (S4) is the 
symmetric matrix whose inverse, when applied to /, reveals the vector of potentials, x. In order to make 
these ideas our own we must work many, many examples. 

2,1,2 Example 

2.1.2.1 Strang Quartet, Step 1 

With respect to the circuit of Figure 2.3 (The fully dressed circuit model), in accordance with step (SI) (list, 
1st bullet, p. 8), we express the six potential differences (always tail minus head) 

ei = x\ — X2 



&2 = X 2 



e-i = x 2 - x 3 



e 4 = x 3 



e5 = X 3 - X4 



e e = x 4 
Such long, tedious lists cry out for matrix representation, to wit e = — (Ax) where 



/ -1 1 



V o 







o \ 






-1 











-1 


1 











-1 











-1 


1 











-1/ 



2.1.2.2 Strang Quartet, Step 2 

Step (S2) (list, 2nd bullet, p. 8), Ohm's Law, states: 

Law 2.1: Ohm's Law 

The current along an edge is equal to the potential drop across the edge divided by the resistance 
of the edge. 

In our case, 



Vj 



Ri 



j = 1,3,5 and yj 



or, in matrix notation, y = Ge where 



G 



R,, 



3 = 2,4,6 



(h 




















1 




















1 
R, 




















1 

Am 




















1 
II, 





I o 














1 

R, 



2.1.2.3 Strang Quartet, Step 3 

Step (S3) (list, 3rd bullet, p. 8), Kirchhoff's Current Law 2 , states: 

Law 2.2: Kirchhoff's Current Law 

The sum of the currents into each node must be zero. 



In our case 



io ~ 2/1 = 



2/1-2/2-2/3 = 



2/3-2/4-2/5 = 



"Kirchhoff's Laws": Section Kirchhoff's Current Law <http://cnx.Org/content/m0015/latest/#current> 
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or, in matrix terms 



where 



B 



2/5 - J/6 = 
By = -f 

( -1 o o o \ 
1-1-10 
1-1-10 

\ 1 -1 ) 



and / 






v o ; 



2.1.2.4 Strang Quartet, Step 4 

Looking back at A: 



( ~ l 


1 





\ 





-1 











-1 


1 











-1 











-1 


1 


V o 








-1/ 



we recognize in B the transpose of A. Calling it such, we recall our main steps 



• (SI) e = -(Ax), 

• (S2) y = Ge, and 



. (S3) A T y=-f. 
On substitution of the first two into the third we arrive, in accordance with (S4) (list, 4th bullet, p. 8), at 

A T GAx = /. (2.1) 

This is a system of four equations for the 4 unknown potentials, x\ through %a- As you know, the system 
(2.1) may have either 1, 0, or infinitely many solutions, depending on / and A T GA. We shall devote (FIX 
ME CNXN TO CHAPTER 3 AND 4) to an unraveling of the previous sentence. For now, we cross our 
fingers and 'solve' by invoking the Matlab program, fibl.m 3 . 



5 http://www.caam.rice.edu/~caam335/cox/lectures/fibl.in 
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Results of a 64 compartment simulation 

Potential along 64 compartment fiber 



T 1 1 - 




Figure 2.4 



Results of a 64 compartment simulation 



Axial current along 64 compartment fiber 



Membrane current along 64 compartment liber 





a ai j a m 04 at os a? &.e oa I 



2 (cm) 



(b) 



Figure 2.5 



This program is a bit more ambitious than the above in that it allows us to specify the number of 
compartments and that rather than just spewing the x and y values it plots them as a function of distance 
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along the fiber. We note that, as expected, everything tapers off with distance from the source and that the 
axial current is significantly greater than the membrane, or leakage, current. 

2,1,3 Example 

We have seen in the previous example (Section 2.1.2: Example) how a current source may produce a potential 
difference across a cell's membrane. We note that, even in the absence of electrical stimuli, there is always a 
difference in potential between the inside and outside of a living cell. In fact, this difference is the biologist's 
definition of 'living.' Life is maintained by the fact that the cell's interior is rich in potassium ions, K + , 
and poor in sodium ions, Na + , while in the exterior medium it is just the opposite. These concentration 
differences beget potential differences under the guise of the Nernst potentials: 

Definition 2.3: Nernst potentials 

RT [Nal RT [KL 

^ = ^ iog kt and E « = -y^w, 

where R is the gas constant, T is temperature, and F is the Faraday constant. 
Associated with these potentials are membrane resistances 

Pm,Na and p m ,K 
that together produce the p m above (list, item 2, p. 5) via 

1 1 1 



Pm Pm,Na Pm,K 

and produce the aforementioned rest potential 



^m Pn 



Pm,Na Pm.Na 



With respect to our old circuit model (Figure 2.3: The fully dressed circuit model), each compartment 
now sports a battery in series with its membrane resistance, as shown in Figure 2.6 (Circuit model with 
resting potentials). 
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Circuit model with resting potentials 




Figure 2.6 



Revisiting steps (Sl-4) (p. 8) we note that in (SI) the even numbered voltage drops are now 

e 2 = x 2 - E m 



e-A = x 3 - E n 



e$ = X4 — E m 
We accommodate such things by generalizing (SI) (list, 1st bullet, p. 8) to: 

• (SI') Express the voltage drops as e = b — Ax where b is the vector of batteries. 
No changes are necessary for (S2) and (S3). The final step now reads, 

• (S4') Combine (SI') (p. 13), (S2) (list, 2nd bullet, p. 8), and (S3) (list, 3rd bullet, p. 
A T GAx = A T Gb + f. 

Returning to Figure 2.6 (Circuit model with resting potentials), we note that 



to produce 



/ 


( \\ 




1 


E m 




1 







\ 


^ 1// 



and A T Gb 



E„ 



( o\ 

1 
1 

\1 J 
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This requires only minor changes to our old code. The new program is called fib2.m 4 and results of its use 
are indicated in the next two figures. 



Results of a 64 compartment simulation with batteries 

Potential along B4 compartment fiber 




Figure 2.7 



4 http://www.caam.rice.edu/~caam335/cox/lectures/fib2.rn 



Results of a 64 compartment simulation with batteries 
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Axial current along e* compartrrent titer 



Membrane current along 64 compartment fiber 





u 11 OS &3 Ofl ;-:. . .= :■ ■■ he 00 I 



z (cm) 



(a) 



(b) 



Figure 2.8 



2.2 CAAM 335 Chapter 1 Exercises 5 

Exercise 2.1 (Solution on p. 17.) 

In order to refresh your matrix- vector multiply skills please calculate, by hand, the product A T GA 
in the 3 compartment case and write out the 4 equations in the vector equation we arrived at in 
step (S4) (2.1): A T GAx = f. 



Exercise 2.2 

We began our discussion with the 'hope' that a multicompartment model could indeed adequately 
capture the fiber's true potential and current profiles. In order to check this one should run fibl.m 6 
with increasing values of N until one can no longer detect changes in the computed potentials. 

• (a) Please run fibl.m 7 with N = 8, 16, 32, and 64. Plot all of the potentials on the same 
(use hold) graph, using different line types for each. (You may wish to alter f ibl .m so that 
it accepts N as an argument). 

Let us now interpret this convergence. The main observation is that the difference equation, (2.2), 
approaches a differential equation. We can see this by noting that 



d(z) 



N 



acts as a spatial 'step' size and that Xk, the potential at (k — 1) d(z), is approximately the value of 
the true potential at (k — 1) d(z). In a slight abuse of notation, we denote the latter 

x((k- l)d(z)) 



5 This content is available online at <http://cnx.Org/content/ml0299/2.8/>. 
6 http:// www.caam.rice.edu/~caam335/cox/lectures/fibl.rn 
7 http:// www.caam.rice.edu/~caam335/cox/lectures/fibl.rn 
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Applying these conventions to (2.2) and recalling the definitions of Ri and R m we see (2.2) become 
7ra 2 -x (0) + 2x (d (z)) - x (2d (z)) 2irad (z) 



Pi d(z 

or, after multiplying through by 



ad(z) ' 

ap m -x (0) + 2x (d [z))-x (2d (z)) 
Pi d(z 2 ) 



-x(d(z)) = 



+ 2x{d{z)) = 



. We note that a similar equation holds at each node (save the ends) and that as N — > oo and 
therefore d (z) — > we arrive at 

d 2 2p. 

-T^X(z) 

dz z ap. 



2 x(z)- -^-x{z)=0 (2.3) 



• (b) With a = ^- show that 

x (z) = asinh ( \[2\±z\ + /3cosh ( \/2p,z) (2.4) 

satisfies (2.3) regardless of a and j3. 
We shall determine a and (3 by paying attention to the ends of the fiber. At the near end we find 

7ra 2 x(0) - x(d(z)) 



Pi d(z) 



'>(} 



which, as d(z) — > becomes 

U^ = - P A (2-5) 

dz Tra z 

At the far end, we interpret the condition that no axial current may leave the last node to mean 

— a; (/) = (2.6) 



• (c) Substitute (2.4) into (2.5) and (2.6) and solve for a and (3 and write out the final x (z). 

• (d) Substitute into x the I, a, pi, and p m values used in fibl.m 8 , plot the resulting function 
(using, e.g., ezplot) and compare this to the plot achieved in part (a) (p. 15). 



i http://www.caam.rice.edu/~caam335/cox/lectures/fibl.rn 
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Solutions to Exercises in Chapter 2 

Solution to Exercise 2.1 (p. 15) 



2.1 Feedback 

The second equation should read 



-X-] + 2X2 — X?, Xo , 

- - (2.7) 



rti ityj 
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Chapter 3 

Matrix Methods for Mechanical Systems 



3.1 A Uniaxial Truss 1 

3.1.1 Introduction 

We now investigate the mechanical prospection of tissue, an application extending techniques developed in 
the electrical analysis of a nerve cell (Section 2.1). In this application, one applies traction to the edges of a 
square sample of planar tissue and seeks to identify, from measurement of the resulting deformation, regions 
of increased 'hardness' or 'stiffness.' For a sketch of the associated apparatus, visit the Biaxial Test site 2 . 

3.1.2 A Uniaxial Truss 



A uniaxial truss 




Figure 3.1 



As a precursor to the biaxial problem (Section 3.2) let us first consider the uniaxial case. We connect 3 
masses with four springs between two immobile walls, apply forces at the masses, and measure the associated 
displacement. More precisely, we suppose that a horizontal force, fj, is applied to each rrij, and produces a 
displacement Xj, with the sign convention that rightward means positive. The bars at the ends of the figure 
indicate rigid supports incapable of movement. The kj denote the respective spring stiffnesses. The analog 



1 This content is available online at <http://cnx.Org/content/ml0146/2.6/>. 
2 http://health. upenn.edu/orl/research/bioengineering/proj2b.jpg 
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of potential difference (see the electrical model (Section 2.1.2.1: Strang Quartet, Step 1)) is here elongation. 
If ej denotes the elongation of the jth spring then naturally, 



ei = xi 



e 2 = x 2 - Xi 



e3 = x 3 - x 2 



e 4 



-X3 



or, in matrix terms, e = Ax, where 



/ 1 \ 

-1 1 

-1 1 

V o o -i y 

We note that ej is positive when the spring is stretched and negative when compressed. This observation, 
Hooke's Law, is the analog of Ohm's Law in the electrical model (Law 2.1, Ohm's Law, p. 9). 

Definition 3.1: Hooke's Law 

1. The restoring force in a spring is proportional to its elongation. We call the constant of propor- 
tionality the stiffness, kj, of the spring, and denote the restoring force by yj. 

2. The mathematical expression of this statement is: yj = kjej, or, 

3. in matrix terms: y = Ke where 

( fci \ 

k 2 

k 3 

fc 4 



K 



The analog of Kirchhoff's Current Law (Section 2.1.2.3: Strang Quartet, Step 3) is here typically called 
'force balance.' 

Definition 3.2: force balance 

1. Equilibrium is synonymous with the fact that the net force acting on each mass must vanish. 

2. In symbols, 

2/i-2/2-/i = 

2/2-2/3-/2 = 



3. or, in matrix terms, By = f where 

h 
\h J 



2/3-2/4-/3 = 



and B 



/ 1 -1 \ 

1-10 
\ 1 -1 / 
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As in the electrical example (Section 2.1.2.4: Strang Quartet, Step 4) we recognize in B the transpose of 
A. Gathering our three important steps: 

e = Ax (3.1) 



y = Ke 



A T y = f 



we arrive, via direct substitution, at an equation for x. Namely, 

(A T y = /) =► (A T Ke = f) => (A T KAx = f) 
Assembling A T KAx we arrive at the final system: 



( ki + k 2 -k 2 

-k 2 k 2 + k 3 -k 3 



\(*A fh\ 



V o 



x 2 



-k 3 k 3 + k A J \ x 3 J 



h 

\h ) 



(3.2) 
(3.3) 



(3.4) 



3.1.3 Gaussian Elimination and the Uniaxial Truss 

Although Matlab solves systems like the one above with ease our aim here is to develop a deeper understand- 
ing of Gaussian Elimination and so we proceed by hand. This aim is motivated by a number of important 
considerations. First, not all linear systems have solutions and even those that do do not necessarily possess 
unique solutions. A careful look at Gaussian Elimination will provide the general framework for not only 
classifying those systems that possess unique solutions but also for providing detailed diagnoses of those 
defective systems that lack solutions or possess too many. 

In Gaussian Elimination one first uses linear combinations of preceding rows to eliminate nonzeros below 
the main diagonal and then solves the resulting triangular system via back-substitution. To firm up our 
understanding let us take up the case where each kj = 1 and so (3.4) takes the form 



/ 2 -1 \ 



-1 2 



(*\ ( h \ 

h 



■T2 



\ -1 2 J \ x 3 J 



We eliminate the (2,1) (row 2, column 1) element by implementing 



(3.5) 



bringing 



new row 2 = old row 2 -\ — row 1 

2 



/ 2 -1 \ / Xl \ ( h \ 



x 2 



V -1 2 J \ x 3 J 
We eliminate the current (3,2) element by implementing 



f2+f 

V h J 



new row 3 = old row 3 H — row 2 
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bringing the upper-triangular system 

( 2 -1 \ / x x \ ( 

-1 X 2 



^ 

u 2 

y o o 



/ 



V*3/ 



2/2 , A 



Kh + 'f + f J 



One now simply reads off 



x 3 



h + 2/2 + 3/3 



This in turn permits the solution of the second equation 

2(x 3 + /2+^) /l+ 2/ 2 + /3 



,r 2 



and, in turn, 



Xi 



3 2 

2:2 + /1 3/i + 2/ 2 + / 3 



2 4 

One must say that Gaussian Elimination has succeeded here. For, regardless of the actual elements of /, we 
have produced an x for which A T KAx = f. 

3.1.4 Alternate Paths to a Solution 

Although Gaussian Elimination remains the most efficient means for solving systems of the form Sx = f it 
pays, at times, to consider alternate means. At the algebraic level, suppose that there exists a matrix that 
'undoes' multiplication by S in the sense that multiplication by 2 _1 undoes multiplication by 2. The matrix 
analog of 2 _1 2 = 1 is 

s-\s = I 

where / denotes the identity matrix (all zeros except the ones on the diagonal). We refer to 6 1-1 as: 

Definition 3.3: Inverse of S 

Also dubbed "S inverse" for short, the value of this matrix stems from watching what happens 
when it is applied to each side of Sx = f. Namely, 

(Sx = /) =* (S^Sx = S- 1 /) =» (Ix = 5- 1 /) => (x = S- 1 /) 



Hence, to solve Sx = f for x it suffices to multiply / by the inverse of S. 

3.1.5 Gauss-Jordan Method: Computing the Inverse of a Matrix 

Let us now consider how one goes about computing S 1-1 . In general this takes a little more than twice the 
work of Gaussian Elimination, for we interpret 

S^S = I 

as n (the size of S) applications of Gaussian elimination, with / running through n columns of the identity 
matrix. The bundling of these n applications into one is known as the Gauss-Jordan method. Let us 
demonstrate it on the S appearing in (3.5). We first augment S with I. 

(2 -1 I 1 \ 

-1 2 -1 I 1 
\ -1 2 J 1 J 
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We then eliminate down, being careful to address each of the three / vectors. This produces 



( 2 









Now, rather than simple back-substitution we instead eliminate up. Eliminating first the (2,3) element we 
find 



(2 

| 
\ 
Now, eliminating the (1,2) element we achieve 





1 







\ 



I i o 



3 
4 

1 / 



1 / 



In the final step we scale each row in order that the matrix on the left takes on the form of the identity. This 
requires that we multiply row 1 by I, row 2 by I, and row 3 by |, with the result 



/ 


1 













i 





V 








1 



1 

2 



Now in this transformation of S into / we have, ipso facto, transformed / to S 1-1 ; i.e., the matrix that 
appears on the right after applying the method of Gauss- Jordan is the inverse of the matrix that began on 
the left. In this case, 

5" x = 



U 






One should check that S 1 / indeed coincides with the x computed above. 

3.1.6 Invertibility 

Not all matrices possess inverses: 

Definition 3.4: singular matrix 

A matrix that does not have an inverse. 
Example 

A simple example is: 

1 1 
1 1 



Alternately, there are 

Definition 3.5: Invertible, or Nonsingular Matrices 

Matrices that do have an inverse. 
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Example 

The matrix S that we just studied is invertible. Another simple example is 



1 

1 1 



3.2 A Small Planar Truss 3 

We return once again to the biaxial testing problem, introduced in the uniaxial truss module (Section 3.1.1: 
Introduction). It turns out that singular matrices are typical in the biaxial testing problem. As our initial 
step into the world of such planar structures let us consider the simple truss in the figure of a simple swing 
(Figure 3.2: A simple swing). 



A simple swing 




Figure 3.2 



We denote by x\ and x 2 the respective horizontal and vertical displacements of m\ (positive is right and 
down) . Similarly, f\ and fi will denote the associated components of force. The corresponding displacements 
and forces at m-i will be denoted by X3, X4 and fs, f±. In computing the elongations of the three springs we 
shall make reference to their unstretched lengths, L\, L 2 , and L 3 . 

Now, if spring 1 connects {0, — Li} to {0,1} when at rest and {0,—Li} to {2:1,3:2} when stretched then 
its elongation is simply 



ei 



^2x 1 2 +{x 2 + L 1 y 



Li 



(3.6) 



3 This content is available online at <http://cnx.Org/content/ml0147/2.6/>. 
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The price one pays for moving to higher dimensions is that lengths are now expressed in terms of square 
roots. The upshot is that the elongations are not linear combinations of the end displacements as they were 
in the uniaxial case (Section 3.1.2: A Uniaxial Truss). If we presume, however, that the loads and stiffnesses 
are matched in the sense that the displacements are small compared with the original lengths, then we may 
effectively ignore the nonlinear contribution in (3.6). In order to make this precise we need only recall the 

Rule 3.1: Taylor development of the square root of (1 + t) 
The Taylor development of \J2\ + t about t = is 

y/21 + t = 1 + - + O (t 2 ) 

where the latter term signifies the remainder. 
With regard to e\ this allows 

ei = \j2x-i 2 + x 2 2 + 2x 2 L x + L 2 - L x 



(3.7) 



(3.1 



If we now assume that 

X\ 2 + X2 2 

is small compared to x 2 (3-9) 

2L\ 

then, as the O term is even smaller, we may neglect all but the first terms in the above and so arrive at 

ei = x 2 
To take a concrete example, if L\ is one meter and x\ and x 2 are each one centimeter, then x 2 is one hundred 

2 , 2 

times Xi 2 l ■ 

With regard to the second spring, arguing as above, its elongation is (approximately) its stretch along 
its initial direction. As its initial direction is horizontal, its elongation is just the difference of the respective 
horizontal end displacements, namely, 

e 2 = x 3 - xi 

Finally, the elongation of the third spring is (approximately) the difference of its respective vertical end 
displacements, i.e., 

e 3 = Xi 

We encode these three elongations in 

( 1 \ 

-10 10 



Ax where 



\ 1/ 



Hooke's law (Definition: "Hooke's Law", p. 20) is an elemental piece of physics and is not perturbed by 
our leap from uniaxial to biaxial structures. The upshot is that the restoring force in each spring is still 
proportional to its elongation, i.e., yj = kjej where kj is the stiffness of the jth spring. In matrix terms, 



y = Ke where K 



( fci \ 

k 2 
\ k 3 J 
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Balancing horizontal and vertical forces at m\ brings 

(-2/2) - /1 = 

and 

2/1-/2 = 

while balancing horizontal and vertical forces at m<i brings 

2/2 - h = 
and 



2/3-/4 = 



We assemble these into 



By = f where B 



( -1 \ 
1 
1 

and recognize, as expected, that B is nothing more than A T . Putting the pieces together, we find that x 
must satisfy Sx = f where 

' ' » ,2 

fci 

-k 2 fc 2 

\ k 3 J 

Applying one step of Gaussian Elimination (Section 3.1.3: Gaussian Elimination and the Uniaxial Truss) 
brings 



(hi -ki \ 



S= A 1 KA 



( k 2 



-ki \ 



fci 


V fc 3 y 



and back substitution delivers 



X2 
X3 

\ X 4 /' 

k 

h 

= /l + /3 



/ 2 
/1 + /3 

V A J 



X4 



X2 



k 

ki 



xi - x 3 



k 

k 2 



The second of these is remarkable in that it contains no components of x. Instead, it provides a condition 
on /. In mechanical terms, it states that there can be no equilibrium unless the horizontal forces on the two 
masses are equal and opposite. Of course one could have observed this directly from the layout of the truss. 
In modern, three-dimensional structures with thousands of members meant to shelter or convey humans 
one should not however be satisfied with the 'visual' integrity of the structure. In particular, one desires 
a detailed description of all loads that can, and, especially, all loads that can not, be equilibrated by the 
proposed truss. In algebraic terms, given a matrix S, one desires a characterization of 
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1. all those / for which Sx = f possesses a solution 

2. all those / for which Sx = f does not possess a solution 

We will eventually provide such a characterization in our later discussion of the column space (Section 4.1) 
of a matrix. 

Supposing now that j\ + / 3 = we note that although the system above is consistent it still fails to 
uniquely determine the four components of x. In particular, it specifies only the difference between x\ and 
X3. As a result both 

/ \ 



k 2 

h 
fei 





V fe/ 



and x 



satisfy Sx = f. In fact, one may add to either an arbitrary multiple of 

/l\ 


1 



(3.10) 



and still have a solution of Sx = f. Searching for the source of this lack of uniqueness we observe some 
redundancies in the columns of S. In particular, the third is simply the opposite of the first. As S is simply 
A T KA we recognize that the original fault lies with A, where again, the first and third columns are opposites. 
These redundancies are encoded in z in the sense that 

Az = 

Interpreting this in mechanical terms, we view z as a displacement and Az as the resulting elongation. In 

Az = 

we see a nonzero displacement producing zero elongation. One says in this case that the truss deforms 
without doing any work and speaks of z as an unstable mode. Again, this mode could have been observed 
by a simple glance at Figure 3.2 (A simple swing). Such is not the case for more complex structures and so 
the engineer seeks a systematic means by which all unstable modes may be identified. We shall see later 
that all these modes are captured by the null space (Section 4.2) of A. 
From 

Sz = 

one easily deduces that S is singular (Definition: "singular matrix", p. 23). More precisely, if S 1-1 were to 
exist then S~ 1 Sz would equal S _1 0, i.e., z = 0, contrary to (3.10). As a result, Matlab will fail to solve 
Sx = f even when / is a force that the truss can equilibrate. One way out is to use the pseudo-inverse, 
as we shall see in the General Planar Truss (Section 3.3) module. 



3.3 The General Planar Truss 4 

Let us now consider something that resembles the mechanical prospection problem introduced in the intro- 
duction to matrix methods for mechanical systems (Section 3.1.1: Introduction). In the figure below we offer 
a crude mechanical model of a planar tissue, say, e.g., an excised sample of the wall of a vein. 



This content is available online at <http://cnx.Org/content/ml0148/2.9/>. 
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A crude tissue model 




Figure 3.3 



Elastic fibers, numbered 1 through 20, meet at nodes, numbered 1 through 9. We limit our observation to 
the motion of the nodes by denoting the horizontal and vertical displacements of node j by x 2 j-i (horizontal) 
and x 2 j (vertical), respectively. Retaining the convention that down and right are positive we note that the 
elongation of fiber 1 is 

ei = x 2 - x 8 



while that of fiber 3 is 



S3 = x 3 - X\. 



As fibers 2 and 4 are neither vertical nor horizontal their elongations, in terms of nodal displacements, are 

not so easy to read off. This is more a nuisance than an obstacle however, for noting our discussion of 

elongation (p. 25) in the small planar truss module, the elongation is approximately just the stretch along 

its undeformed axis. With respect to fiber 2, as it makes the angle — J with respect to the positive horizontal 

axis, we find 

7T \ . / n\ x 9 - xi + x 2 - x 10 

I -.TinSlTl I I = 



62 = SgCOS 



cciosin 



Similarly, as fiber 4 makes the angle — ^ with respect to the positive horizontal axis, its elongation is 



3tt 

"T 



e4 = X7COS I — — I — xgsin 
These are both direct applications of the general formula 

e,- = X2n-lCOS(9j) - 



x 3 - x 7 + x 4 - x% 



22 



x 2n sm (tlj) 



(3.11) 



for fiber j, as depicted in Figure 3.4 below, connecting node m to node n and making the angle 9j with 
the positive horizontal axis when node m is assumed to lie at the point (0,0). The reader should check that 
our expressions for e\ and e% indeed conform to this general formula and that e 2 and e\ agree with ones 
intuition. For example, visual inspection of the specimen suggests that fiber 2 can not be supposed to stretch 
(i.e., have positive e 2 ) unless xg > x\ and/or x 2 > xiq. Does this jive with (3.11)? 
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original 




deformed 



2m-l 



Figure 3.4: Elongation of a generic bar, see (3.11) 



Applying (3.11) to each of the remaining fibers we arrive at e = Ax where A is 20-by-18, one row for 
each fiber, and one column for each degree of freedom. For systems of such size with such a well defined 
structure one naturally hopes to automate the construction. We have done just that in the accompanying 
M-file 5 and diary 6 . The M-file begins with a matrix of raw data that anyone with a protractor could have 
keyed in directly from Figure 3.3 (A crude tissue model): 



data = 

1 
1 
1 
2 



-pi/2 
-pi/4 
) 
-3*pi/4 







. . . and so on. 



7, one row of data for each fiber, the 
'/, first two columns are starting and ending 
7, node numbers, respectively, while the third is the 
7, angle the fiber makes with the positive horizontal axis 



] 



This data is precisely what (3.11) requires in order to know which columns of A receive the proper cos or 
sin. The final A matrix is displayed in the diary 7 . 

The next two steps are now familiar. If K denotes the diagonal matrix of fiber stiffnesses and / denotes 
the vector of nodal forces then 

y = Ke and A T y = f 

and so one must solve Sx = f where S = A T KA. In this case there is an entire three-dimensional class of 
z for which Az = and therefore Sz = 0. The three indicates that there are three independent unstable 
modes of the specimen, e.g., two translations and a rotation. As a result S is singular and x = S\f in 
MATLAB will get us nowhere. The way out is to recognize that S has 18 — 3 = 15 stable modes and that 
if we restrict S to 'act' only in these directions then it 'should' be invertible. We will begin to make these 
notions precise in discussions on the Fundamental Theorem of Linear Algebra. (Section 4.5) For now let us 



5 http 
6 http 
7 http 



//cnx.rice.edu/modules/ml0148/latest /fiber. m 
//cnx.rice.edu/modules/ml0148/latest/lec2adj 
//cnx.rice.edu/modules/ml0148/latest/lec2adj 
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note that every matrix possesses such a pseudo-inverse and that it may be computed in MATLAB via the 
pinv command. Supposing the fiber stiffnesses to each be one and the edge traction to be of the form 



/=( -1 1 1 1 1 -1 1 -1 -1 -1 1 -1 ) 

we arrive at x via x=pinv(S)*f and offer below its graphical representation. 

3.3.1 Before- After Plot 



T 







o o 






Figure 3.5: Before and after shots of the truss in Figure 3.3 (A crude tissue model). The solid (dashed) 
circles correspond to the nodal positions before (after) the application of the traction force, /. 



3.4 CAAM 335 Chapter 2 Exercises 8 

Exercise 3.1 

With regard to the uniaxial truss figure (Figure 3.1: A uniaxial truss), 

• (i) Derive the A and K matrices resulting from the removal of the fourth spring, 

• (ii) Compute the inverse, by hand via Gauss- Jordan (Section 3.1.5: Gauss- Jordan Method: 
Computing the Inverse of a Matrix), of the resulting A T KA with fci = &2 = &3 = k 

• (iii) Use the result of (ii) to find the displacement corresponding to the load / = (0, 0, F) . 

Exercise 3.2 

Generalize example 3, the general planar truss (Section 3.3), to the case of 16 nodes connected by 
42 fibers. Introduce one stiff (say k = 100) fiber and show how to detect it by 'properly' choosing /. 
Submit your well-documented M-file as well as the plots, similar to the before-after plot (Figure 3.6) 
in the general planar module (Section 3.3), from which you conclude the presence of a stiff fiber. 



3 This content is available online at <http://cnx.Org/content/ml0300/2.6/>. 
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Figure 3.6: A copy of the before-after figure from the general planar module. 
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Chapter 4 

The Fundamental Subspaces 



4.1 Column Space 1 
4.1.1 The Column Space 

We begin with the simple geometric interpretation of matrix-vector multiplication. Namely, the multipli- 
cation of the n-by-1 vector x by the m-by-n matrix A produces a linear combination of the columns of A. 
More precisely, if aj denotes the jth column of A, then 



Ax 



( ai 



d-2 



( -1 \ 

X 2 



(4.1) 



\ x n J 



x\a\ + x 2 a 2 H + x n a r , 



The picture that I wish to place in your mind's eye is that Ax lies in the subspace spanned (Definition: 
"Span", p. 44) by the columns of A. This subspace occurs so frequently that we find it useful to distinguish 
it with a definition. 

Definition 4.1: Column Space 

The column space of the m-by-n matrix S is simply the span of the its columns, i.e. Ra (S) = 
{Sx | x € R n }. This is a subspace (Section 4.7.2) of 5R m . The notation Ra stands for range in this 
context. 



4.1.2 Example 

Let us examine the matrix: 



/ 1 \ 

-10 10 
\ 1/ 



(4.2) 



lr This content is available online at <http://cnx.Org/content/ml0266/2.9/>. 
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The column space of this matrix is: 

' / \ /i\ /o\ /o\ 

Ra(A) =lxi -1 + x 2 + x 3 1 + x 4 
As the third column is simply a multiple of the first, we may write: 



Ra(A) = < 



' 


(°) 




(l) 




(°) 


Xi 


1 


+ x 2 





+ x 3 







V ° / 




V°/ 




V i/ 



X £ 



(4.3) 



X £ 



(4.4) 



As the three remaining columns are linearly independent (Definition: "Linear Independence", p. 44) we 
may go no further. In this case, Ra (A) comprises all of R 3 . 

4,1,3 Method for Finding a Basis 

To determine the basis for Ra (A) (where A is an arbitrary matrix) we must find a way to discard its 
dependent columns. In the example above, it was easy to see that columns 1 and 3 were colinear. We seek, 
of course, a more systematic means of uncovering these, and perhaps other less obvious, dependencies. Such 
dependencies are more easily discerned from the row reduced form. In the reduction of the above problem, 
we come very easily to the matrix 

( -1 1 \ 



10 
\ 1 J 



(4.5) 



Once we have done this, we can recognize that the pivot column are the linearly independent columns of 
A re( j. One now asks how this might help us distinguish the independent columns of A. For, although the 
rows of j4 rc d are linear combinations of the rows of A, no such thing is true with respect to the columns. 
The answer is: pay attention to the indices of the pivot columns. In our example, columns {1, 2, 4} 
are the pivot columns of A rc( j and hence the first, second, and fourth columns of A, i.e., 



(4.6) 



\( ° ) 




(i) 




(°) 


-1 


i 





t 





IA ° ) 




\o) 




V i/ 



comprise a basis for Ra (A). In general: 

Definition 4.2: A Basis for the Column Space 



Suppose A is m-by-n. If columns { < 



l,...,r} are the pivot columns of A re( j then columns 



{cj | j = 1, ..., r } of A constitute a basis for Ra (A). 



4.2 Null Space 2 
4.2.1 Null Space 

Definition 4.3: Null Space 

The null space of an m-by-n matrix A is the collection of those vectors in 



that A maps to the 



2 This content is available online at <http://cnx.Org/content/ml0293/2.9/>. 



zero vector in W 71 . More precisely, 



JV(A) = {xeR n \Ax = 0} 
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4.2.2 Null Space Example 

As an example, we examine the matrix A: 



( 1 \ 

-10 10 
\ 1/ 



It is fairly easy to see that the null space of this matrix is: 

f /l\ 



,yV{A)={ 




1 



te 



(4.7) 



(41 



This is a line in M 4 . 

The null space answers the question of uniqueness of solutions to Sx = /. For, if Sx = / and Sy = f 
then S (x — y) = Sx — Sy = f — f = and so (x — y) € JV (S). Hence, a solution to Sx = f will be unique 
if, and only if, ^V (S) = {0}. 

4.2.3 Method for Finding the Basis 

Let us now exhibit a basis for the null space of an arbitrary matrix A. We note that to solve Ax = is to 
solve A Te( ±x = 0. With respect to the latter, we suppose that 



{ Cj \j = {!,.. .,r}} 

are the indices of the pivot columns and that 

{ Cj \j = {r+l,...,n}} 

are the indices of the nonpivot columns. We accordingly define the r pivot variables 

{x Cj \j = {!,..., r}} 

and the n — r free variables 

{x Cj \ j = {r+l,...,n}} 



(4.9) 
(4.10) 

(4.11) 

(4.12) 



One solves A ro ax = by expressing each of the pivot variables in terms of the nonpivot, or free, variables. 
In the example above, xi, xi, and x\ are pivot while xj, is free. Solving for the pivot in terms of the free, we 
find Xi = 0, X3 = x\, xi = 0, or, written as a vector, 



x = x 3 




1 



(4.13) 
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where X3 is free. As X3 ranges over all real numbers the x above traces out a line in R 4 . This line is precisely 
the null space of A. Abstracting these calculations we arrive at: 

Definition 4.4: A Basis for the Null Space 

Suppose that A is m-by-n with pivot indices {cj \ j = {1, . . . , r}} and free indices 
{cj I j = {r + 1, . . . , n} }. A basis for jY {A) may be constructed of n — r vectors {z 1 , z 2 , . . . , z n ~ r } 
where z k , and only z k , possesses a nonzero in its c r+ k component. 

4.2.4 A MATLAB Observation 

As usual, MATLAB has a way to make our lives simpler. If you have defined a matrix A and want to 
find a basis for its null space, simply call the function null (A) . One small note about this function: if one 
adds an extra flag, 'r', as in null (A, 'r'), then the basis is displayed "rationally" as opposed to purely 
mathematically. The MATLAB help pages define the difference between the two modes as the rational mode 
being useful pedagogically and the mathematical mode of more value (gasp!) mathematically. 

4.2.5 Final thoughts on null spaces 

There is a great deal more to finding null spaces; enough, in fact, to warrant another module. One important 
aspect and use of null spaces is their ability to inform us about the uniqueness of solutions. If we use the 
column space 3 to determine the existence of a solution x to the equation Ax = b. Once we know that a 
solution exists it is a perfectly reasonable question to want to know whether or not this solution is the only 
solution to this problem. The hard and fast rule is that a solution x is unique if and only if the null space of 
A is empty. One way to think about this is to consider that if Ax = does not have a unique solution then, 
by linearity, neither does Ax = b. Conversely, if (Az = 0) and (z 7^ 0) and (Ay = b) then A (z + y) = b 
as well. 

4.3 The Null and Column Spaces: An Example 4 

4.3.1 Preliminary Information 

Let us compute bases for the null and column spaces of the adjacency matrix associated with the ladder 
below. 

3 <http://cnx.org/content/columnspace/latest/> 

4 This content is available online at <http://cnx.Org/content/ml0368/2.4/>. 
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An unstable ladder? 



2 

■O 



4 



8 



Figure 4.1 



The ladder has 8 bars and 4 nodes, so 8 degrees of freedom. Denoting the horizontal and vertical 
displacements of node j by X2j-i and X2j, respectively, we arrive at the A matrix 



( 1 
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-1 





1 
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-1 











1 

















-1 
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-1 
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-1 


0/ 



4,3.2 Finding a Basis for the Column Space 

To determine a basis for 1Z (A) we must find a way to discard its dependent columns. A moment's reflection 
reveals that columns 2 and 6 are colinear, as are columns 4 and 8. We seek, of course, a more systematic 
means of uncovering these and perhaps other less obvious dependencies. Such dependencies are more easily 
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discerned from the row reduced form (Section 4.7.3) 



*rcd 



rref (A) 



1 1 




















\ 





1 











-1 














1 


























1 











-1 














1 
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v° 




















o ) 



Recall that rref performs the elementary row operations necessary to eliminate all nonzeros below the 
diagonal. For those who can't stand to miss any of the action I recommend rref movie. 

Each nonzero row of A xc & is called a pivot row. The first nonzero in each row of A re( j is called a pivot. 
Each column that contains a pivot is called a pivot column. On account of the staircase nature of A re ^ we 
find that there are as many pivot columns as there are pivot rows. In our example there are six of each and, 
again on account of the staircase nature, the pivot columns are the linearly independent (Definition: "Linear 
Independence", p. 44) columns of A rc( ±. One now asks how this might help us distinguish the independent 
columns of A. For, although the rows of A re d are linear combinations of the rows of A, no such thing is true 
with respect to the columns. In our example, columns {1,2,3,4,5,7} are the pivot columns. In general: 

Proposition 4.1: 

Suppose A is m-by-n. If columns {cj \ j = 1, ...,r} are the pivot columns of j4 rcc j. then columns 
{cj \ j = 1, ..., r } of A constitute a basis for 1Z (A). 
Proof: 

Note that the pivot columns of A re( j. are, by construction, linearly independent. Suppose, however, 
that columns { Cj \ j = l,...,r} of A are linearly dependent. In this case there exists a nonzero 
x e W 1 for which Ax = and 



Xk 



, k<£{cj\j= l,...,r} 



(4.14) 



Now Ax = necessarily implies that ^cd^ = 0, contrary to the fact that columns { c } ■. \ j = 1, ..., r } 
are the pivot columns of A vc d- 

We now show that the span of columns {cj \ j = l,...,r} of A indeed coincides with 1Z (A). 
This is obvious if r = n, i.e., if all of the columns are linearly independent. If r < n, there exists a 
1 & i c j I 3 ' = 1; •••7 r }- Looking back at A re( j. we note that its qth column is a linear combination 
of the pivot columns with indices not exceeding q. Hence, there exists an x satisfying (4.14) and 
A ro ax = 0, and x q = 1. This x then necessarily satisfies Ax = 0. This states that the qth column 
of A is a linear combination of columns {cj \ j = 1, ...,r} of A. 



4,3.3 Finding a Basis for the Null Space 

Let us now exhibit a basis for JY (A). We exploit the already mentioned fact that JV (A) 
Regarding the latter, we partition the elements of x into so called pivot variables, 



^(A-ed)- 



{x Cj | j = l,-,r} 



and free variables 



{x k | k <£ {cj | j = l,...,r}} 



39 

There are evidently n — r free variables. For convenience, let us denote these in the future by 

{x Cj \j = r+ l,...,n} 

One solves A YC( ±x = 0, by expressing each of the pivot variables in terms of the nonpivot, or free, variables. 
In the example above, x\, X2, 23, 24, 25, and 27 are pivot while xq and x% are free. Solving for the pivot in 
terms of the free we find 

27 = 



x 5 = 



X4 = x$ 



23 = 



2 2 = 2 6 



or, written as a vector, 



2 = 26 



21 =0 



/o\ /o\ 



2 8 



\o/ 



V 1/ 



(4.15) 



where x§ and 2s are free. As xq and 2s range over all real numbers, the 2 above traces out a plane in M 8 . 
This plane is precisely the null space of A and (4.15) describes a generic element as the linear combination 
of two basis vectors. Compare this to what MATLAB returns when faced with null (A, 'r')- Abstracting 
these calculations we arrive at 



Proposition 4.2: 

Suppose that A is m-by-n with pivot indices {cj \ j = 1, ...,r} and 
{cj j j = r + 1, ...,n}. A basis for Jf (A), may be constructed of n — r vectors [z 
where z k , and only z k , possesses a nonzero in its c r+ k component. 



free 



indices 

;Z n ~ r } 



4,3.4 The Physical Meaning of Our Calculations 

Let us not end on an abstract note however. We ask what H (A) and Ji (A) tell us about the ladder. 
Regarding H (A) the answer will come in the next chapter. The null space calculation however has revealed 
two independent motions against which the ladder does no work! Do you see that the two vectors in (4.15) 
encode rigid vertical motions of bars 4 and 5 respectively? As each of these lies in the null space of A, the 
associated elongation is zero. Can you square this with the ladder as pictured in Figure 4.1 (An unstable 
ladder?)? I hope not, for vertical motion of bar 4 must 'stretch' bars 1, 2, 6, and 7. How does one resolve 
this (apparent) contradiction? 
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4.4 Left Null Space 5 
4.4.1 Left Null Space 

If one understands the concept of a null space (Section 4.2), the left null space is extremely easy to understand. 

Definition 4.5: Left Null Space 

The Left Null Space of a matrix is the null space (Section 4.2) of its transpose, i.e., 

rf (A T ) = { y e R m | A T y = } 

The word "left" in this context stems from the fact that A T y = is equivalent to y T A = where y 
"acts" on A from the left. 



4.4.2 Example 

:ic 
null space of A T . If 



As A re d was the key to identifying the null space (Section 4.2) of A, we shall see that Aj cd is the key to the 



then 



and so 



red 



(l l\ 

1 2 
\1 3/ 



1 1 1 

1 2 3 

1 1 1 

1 2 



(4.16) 



(4.17) 



(4.18) 



We solve ^4^, d = by recognizing that y\ and yi are pivot variables while 2/3 is free. Solving Aj cd y = for 
the pivot in terms of the free we find yi = — (22/3) and y\ = 2/3 hence 





' 


( 1 \ 


jY (A T ) = < 


2/3 


-2 






V 1 / 



2/3 e 



(4.19) 



4.4.3 Finding a Basis for the Left Null Space 

The procedure is no different than that used to compute the null space (Section 4.2) of A itself. In fact 

Definition 4.6: A Basis for the Left Null Space 

Suppose that A T is n-by-m with pivot indices {cj \ j = {1, . . . ,r}} and free indices 
{cj \ j = {r + 1, . . . ,m}}. A basis for jV (A T ) may be constructed of m — r vectors 



{z\. 



r } where z k , and only z k , possesses a nonzero in its c r +k component. 



5 This content is available online at <http://cnx.Org/content/ml0292/2.7/>. 
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4.5 Row Space 6 
4,5.1 The Row Space 

As the columns of A T are simply the rows of A we call Ra (A T ) the row space of A T . More precisely 

Definition 4.7: Row Space 

The row space of the m-by-n matrix A is simply the span of its rows, i.e., 

Ra (A T ) = { A T y | y e R m } 
This is a subspace of R n . 



4.5.2 Example 

Let us examine the matrix: 



The row space of this matrix is: 



ftH (A T ) = < 



( 1 \ 

-10 10 
\ 1/ 



' 


f°) 




( -l) 




f °\ 


2/1 


1 



+ 2/2 



1 


+ 2/3 








{° J 




^ o J 




^ 1/ 



i/e 



(4.20) 



(4.21) 



As these three rows are linearly independent (Definition: "Linear Independence", p. 44) we may go no 
further. We "recognize" then 1Z~\ (A T ) as a three dimensional subspace (Section 4.7.2) of K 4 . 

4.5.3 Method for Finding the Basis of the Row Space 

Regarding a basis for "RA (A T ) we recall that the rows of A le d, the row reduced form (Section 4.7.3) of the 
matrix A, are merely linear combinations of the rows of A and hence 



ftH (A T ) = ftH (A Icd ) 

This leads immediately to: 

Definition 4.8: A Basis for the Row Space 

Suppose A is m-by-n. The pivot rows of A re( j constitute a basis for 1Z~\ (^4 T ). 
With respect to our example, 

r/ \ /-i\ /oM 



i 







1 

\ o / 






\ 1/J 



comprises a basis for 1Z~\ (A T ). 



(4.22) 



(4.23) 



6 This content is available online at <http://cnx.Org/content/ml0296/2.7/>. 
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4.6 Exercises: Columns and Null Spaces 7 



Exercises 

1. I encourage you to use rref and null for the following. 

• (i) Add a diagonal crossbar between nodes 3 and 2 in the unstable ladder figure (Figure 4.1: An 
unstable ladder?) and compute bases for the column and null spaces of the new adjacency matrix. 
As this crossbar fails to stabilize the ladder, we shall add one more bar. 

• (ii) To the 9 bar ladder of (i) add a diagonal cross bar between nodes 1 and the left end of bar 6. 
Compute bases for the column and null spaces of the new adjacency matrix. 

2. We wish to show that N (A) = N (A T A) regardless of A. 

• (i) We first take a concrete example. Report the findings of null when applied to A and A T A 
for the A matrix associated with the unstable ladder figure (Figure 4.1: An unstable ladder?). 

• (ii) Show that N (A) C N (A T A), i.e. that if Ax = then A T Ax = 0. 

• (iii) Show that N (A T A) C N {A), i.e., that if A T Ax = then Ax = 0. (Hint: if A T Ax = then 
x T A T Ax = 0.) 



3. Suppose that A is m-by-n and that N (A) 



Argue that A must be the zero matrix. 



4.7 Appendices/Supplements 



4.7.1 Vector Space 8 

4.7.1.1 Introduction 

You have long taken for granted the fact that the set of real numbers, R, is closed under addition and 
multiplication, that each number has a unique additive inverse, and that the commutative, associative, and 
distributive laws were right as rain. The set, C, of complex numbers also enjoys each of these properties, as 
do the sets R n and C" of columns of n real and complex numbers, respectively. 

To be more precise, we write x and y in WL n as 

x= (x 1 ,x 2 ,.-.,x n ) T 

y= {yi,V2,---,Vn) 

and define their vector sum as the elementwise sum 



x + y 



( xi + yi \ 

X2 + 2/2 



(4.24) 



\ x n + y n J 
and similarly, the product of a complex scalar, z € C with x as: 

/ zxi \ 

ZX 2 

\ zx n J 



(4.25) 



7 This content is available online at <http://cnx.Org/content/ml0367/2.4/>. 
8 This content is available online at <http://cnx.Org/content/ml0298/2.6/>. 
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4.7.1.2 Vector Space 

These notions lead naturally to the concept of vector space. A set V is said to be a vector space if 



x for each x and y in V 

x + y + z for each x, y, and z in V 



1. x + y = y 

2. x + y + z - 

3. There is a unique "zero vector" such that x + = x for each linV 

4. For each x in V there is a unique vector —x such that x H — x = 0. 

5. la; = x 

6. (C1C2) x = c\ (cix) for each xh\V and ci and C2 in C. 

7. c (a; + y) = ex + cy for each x and y in V and c in C. 

8. (ci + C2) a; = C\X + Cix for each x in V and ci and C2 in C. 



4,7.2 Subspaces 9 

4.7.2.1 Subspace 

A subspace is a subset of a vector space (Section 4.7.1) that is itself a vector space. The simplest example 
is a line through the origin in the plane. For the line is definitely a subset and if we add any two vectors on 
the line we remain on the line and if we multiply any vector on the line by a scalar we remain on the line. 
The same could be said for a line or plane through the origin in 3 space. As we shall be travelling in spaces 
with many many dimensions it pays to have a general definition. 

Definition 4.9: A subset S of a vector space V is a subspace of V when 

1. if x and y belong to S then so does x + y. 

2. if x belongs to S and t is real then tx belong to S. 

As these are oftentimes unwieldy objects it pays to look for a handful of vectors from which the entire 
subset may be generated. For example, the set of x for which x\ + x-i + X3 + 24 = constitutes a subspace 
of M 4 . Can you 'see' this set? Do you 'see' that 

/-l\ 
-1 



V ) 

( -i\ 


1 

V ) 

and 



and 



1 _1 \ 








V 1 / 



not only belong to a set but in fact generate all possible elements? More precisely, we say that these vectors 
span the subspace of all possible solutions. 



9 This content is available online at <http://cnx.Org/content/ml0297/2.6/>. 
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Definition 4.10: Span 

A finite collection {s\, S2, ■ ■ ■ , s n } of vectors in the subspace S is said to span S if each element of 
S can be written as a linear combination of these vectors. That is, if for each s € S there exist n 
reals {x±,X2, ■ • • , x n } such that s = x\S\ + x^s-i + • • • + x n s n . 

When attempting to generate a subspace as the span of a handful of vectors it is natural to ask what is 
the fewest number possible. The notion of linear independence helps us clarify this issue. 

Definition 4.11: Linear Independence 

A finite collection {si, s 2 , ■ ■ ■ , s n } of vectors is said to be linearly independent when the only reals, 
{x\, x-ii • • • , x n } for which x\ + x 2 + • • ■ + x n = are x\ = xi = • • • = x n = 0. In other words, when 
the null space (Section 4.2) of the matrix whose columns are {si, S2, ■ ■ ■ , s n } contains only the zero 
vector. 

Combining these definitions, we arrive at the precise notion of a 'generating set.' 

Definition 4.12: Basis 

Any linearly independent spanning set of a subspace S is called a basis of S. 

Though a subspace may have many bases they all have one thing in common: 

Definition 4.13: Dimension 

The dimension of a subspace is the number of elements in its basis. 



4,7.3 Row Reduced Form 10 

4.7.3.1 Row Reduction 

A central goal of science and engineering is to reduce the complexity of a model without sacrificing its 
integrity. Applied to matrices, this goal suggests that we attempt to eliminate nonzero elements and so 
'uncouple' the rows. In order to retain its integrity the elimination must obey two simple rules. 

Definition 4.14: Elementary Row Operations 

1. You may swap any two rows. 

2. You may add to a row a constant multiple of another row. 

With these two elementary operations one can systematically eliminate all nonzeros below the diagonal. 
For example, given 

1 \ 



/ 




-1 





V i 



1 

2 3 



(4.26) 



it seems wise to swap the first and fourth rows and so arrive at 



/ 



V o 



2 3 4 \ 

1 

1 

1 J 



(4.27) 



°This content is available online at <http://cnx.Org/content/ml0295/2.6/>. 
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adding the first row to the third now produces 

/ 1 2 3 4 \ 

10 

2 4 4 

\ 1 ) 

subtracting twice the second row from the third yields 

/ 1 2 3 4 \ 

10 

4 4 

\ 1 ) 



(4.28) 



(4.29) 



a matrix with zeros below its diagonal. This procedure is not restricted to square matrices. For example, 
given 



/l 1 1 l\ 

2 4 4 2 
\ 3 5 5 3 / 



(4.30) 



we start at the bottom left then move up and right. Namely, we subtract 3 times the first row from the 
third and arrive at 

(l 1 1 



2 4 4 
\ 2 2 
and then subtract twice the first row from the second, 



2 



(4.31) 





1 l l 


1 


1 \ 




2 


2 





\ 2 


2 


0/ 


and finally subtract the second row from the third, 






/l 1 


1 


l\ 




2 


2 







^00 





0/ 



(4.32) 



(4.33) 



It helps to label the before and after matrices. 

Definition 4.15: The Row Reduced Form 

Given the matrix A we apply elementary row operations until each nonzero below the diagonal is 
eliminated. We refer to the resulting matrix as A rc( \. 
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4.7.3.2 Uniqueness and Pivots 

As there is a certain amount of flexibility in how one carries out the reduction it must be admitted that 
the reduced form is not unique. That is, two people may begin with the same matrix yet arrive at different 
reduced forms. The differences however are minor, for both will have the same number of nonzero rows and 
the nonzeros along the diagonal will follow the same pattern. We capture this pattern with the following 
suite of definitions, 

Definition 4.16: Pivot Row 

Each nonzero row of A rc( j is called a pivot row. 

Definition 4.17: Pivot 

The first nonzero term in each row of A Te d is called a pivot. 

Definition 4.18: Pivot Column 

Each column of A ro( j that contains a pivot is called a pivot column. 

Definition 4.19: Rank 

The number of pivots in a matrix is called the rank of that matrix. 

Regarding our example matrices, the first (4.26) has rank 4 and the second (4.30) has rank 2. 

4.7.3.3 Row Reduction in MATLAB 

MATLAB's rref command goes full-tilt and attempts to eliminate ALL off diagonal terms and to leave 
nothing but ones on the diagonal. I recommend you try it on our two examples. You can watch its individual 
decisions by using rref movie instead. 



Chapter 5 

Least Squares 

5.1 Least Squares 1 

5.1.1 Introduction 

We learned in the previous chapter that Ax = b need not possess a solution when the number of rows of A 
exceeds its rank, i.e., r < m. As this situation arises quite often in practice, typically in the guise of 'more 
equations than unknowns,' we establish a rationale for the absurdity Ax = b. 

5.1.2 The Normal Equations 

The goal is to choose x such that Ax is as close as possible to b. Measuring closeness in terms of the sum of 
the squares of the components we arrive at the 'least squares' problem of minimizing 

res 

{\\ Ax - b \\f = {Ax - bf {Ax - b) (5.1) 

over all x £ R. The path to the solution is illuminated by the Fundamental Theorem. More precisely, we 
write b=b R + b N , b R £ R {A) and b N £ N (A T ) . On noting that (i) {Ax - b R ) £ R {A) , x £ R n 
and (ii)M(A) ± N (A T ) we arrive at the Pythagorean Theorem. 

Pythagorean Theorem 

norm 2 {Ax - b) = (|| Ax - b R + b N \\f 

= (HAr-Mlf + dlMI) 2 

It is now clear from the Pythagorean Theorem (5.2: Pythagorean Theorem) that the best x is the one that 
satisfies 

Ax = b R (5.3) 

As b R £ R {A) this equation indeed possesses a solution. We have yet however to specify how one computes 
b R given 6. Although an explicit expression for b R , the so called orthogonal projection of b onto R(^4), 
in terms of A and b is within our grasp we shall, strictly speaking, not require it. To see this, let us note 
that if x satisfies the above equation (5.3) then 

Ax — b = Ax — bo + b\r 

(5.4) 

= -bN 



lr This content is available online at <http://cnx.Org/content/ml0371/2.9/>. 
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As &at is no more easily computed than bn you may claim that we are just going in circles. The 'practical' 
information in the above equation (5.4) however is that {Ax — b) e A T , i.e., A T (Ax — b) = 0, i.e., 



A T Ax = A T b 



(5.5) 



As A T b e K (^4 T ) regardless of b this system, often referred to as the normal equations, indeed has 
a solution. This solution is unique so long as the columns of A T A are linearly independent, i.e., so long 
as N (A T A) = {0}. Recalling Chapter 2, Exercise 2, we note that this is equivalent to N (A) = {0}. We 
summarize our findings in 

Theorem 5.1: 

The set of x G b^ for which the misfit (|| Ax — b \\) is smallest is composed of those x for which 
A T Ax = A T b. There is always at least one such x. There is exactly one such x if N (A) = {0}. 



As a concrete example, suppose with reference to Figure 5.1 that A 



(l l\ 

1 



and b 



1 

V 1/ 




2 -2 



'2 



Figure 5.1: The decomposition of b. 
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As b / K (A) there is no a; such that Ax = b. Indeed, (|| Ax - b \\f = (cci + x 2 + -lf + {x 2 — 1) 2 + 1 > 1, 

I , . 

with the minimum uniquely attained at x = I | , in agreement with the unique solution of the above 

equation (5.5), for A T A = and A T b = \ . We now recognize, a posteriori, that bn = Ax 




1 

v o y 



is the orthogonal projection of 6 onto the column space of A. 



5.1.3 Applying Least Squares to the Biaxial Test Problem 

We shall formulate the identification of the 20 fiber stiffnesses in this previous figure (Figure 3.3: A crude 
tissue model), as a least squares problem. We envision loading, /, the 9 nodes and measuring the associated 
18 displacements, x. From knowledge of x and / we wish to infer the components of K = diag (k) where k 
is the vector of unknown fiber stiffnesses. The first step is to recognize that 

A T KAx = f (5.6) 

may be written as 

Bk = f , B = A T diag {Ax) (5.7) 

Though conceptually simple this is not of great use in practice, for B is 18-by-20 and hence the above 
equation (5.7) possesses many solutions. The way out is to compute k as the result of more than one 
experiment. We shall see that, for our small sample, 2 experiments will suffice. To be precise, we suppose 
that x 1 is the displacement produced by loading f 1 while x 2 is the displacement produced by loading f 2 . 

I A T diag (Ax 1 ) \ I f 1 , 

We then piggyback the associated pieces in B = and / = This B is 36-by-20 

\ A T 'diag (Ax 2 ) j \f 2 

and so the system Bk = f is overdetermined and hence ripe for least squares. 

We proceed then to assemble B and /. We suppose f 1 and f 2 to correspond to horizontal and vertical 
stretching 

f 1 =(-100010-100010-100010) (5.8) 

f 2 =(0101010000000 -10 -lO-l) (5.9) 

respectively. For the purpose of our example we suppose that each kj = 1 except k% = 5. We assemble 
A T KA as in Chapter 2 and solve 

A T KAx 3 = f (5.10) 

with the help of the pseudoinverse. In order to impart some 'reality' to this problem we taint each a;- 7 ' with 
10 percent noise prior to constructing B. Regarding 

B T Bk = B T f (5.11) 

we note that Matlab solves this system when presented with k=B\f when B is rectangular. We have plotted 
the results of this procedure in the Figure 5.2. The stiff fiber is readily identified. 
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10 15 

fiber number 



25 



Figure 5.2: Results of a successful biaxial test. 



5.1.4 Projections 

From an algebraic point of view (5.5) )is an elegant reformulation of the least squares problem. Though easy 
to remember it unfortunately obscures the geometric content, suggested by the word 'projection,' of (5.4). 
As projections arise frequently in many applications we pause here to develop them more carefully. With 
respect to the normal equations we note that if N (A) = {0} then 



{A T A)' 1 A T b 



and so the orthogonal projection of b onto R (A) is: 



b R = Ax 

= A(A T Ay 1 A 1 



Defining 



1 aT 



P = A(A T A) A 



(5.12) 

(5.13) 
(5.14) 
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(5.13) takes the form bn = Pb. Commensurate with our notion of what a 'projection' should be we 
expect that P map vectors not in R (^4) onto M. (A) while leaving vectors already in R (A) unscathed. More 
succinctly, we expect that Pb^ = 6^, i.e., Pb^ = Pb^. As the latter should hold for all 6 € R m we expect 
that 

P 2 = P (5.15) 

With respect to (5.14) we find that indeed 

P 2 = A(A T A)~ 1 A T A(A T A)~ 1 A T 

= A(A T A)' 1 A T (5.16) 

= P 

We also note that the P in (5.14) is symmetric. We dignify these properties through 

Definition 5.1: orthogonal projection 

A matrix P that satisfies P 2 = P is called a projection. A symmetric projection is called an 
orthogonal projection. 

We have taken some pains to motivate the use of the word 'projection.' You may be wondering however 
what symmetry has to do with orthogonality. We explain this in terms of the tautology 

b=Pb-Ib (5.17) 

Now, if P is a projection then so too is J — P. Moreover, if P is symmetric then the dot product of b's two 
constituents is 

(Pbf{I-P)b = b T P T {I-P)b 

= b T (P-P 2 )b 

T V ' (5-18) 

= b T 0b 

= 

i.e., Pb is orthogonal to (/ — P) b. As examples of a nonorthogonal projections we offer P = 
( 1 \ 

and I — P. Finally, let us note that the central formula, P = ^4(^4 T ^4) = A T , 

\^ ? J 
even a bit more general than advertised. It has been billed as the orthogonal projection onto the column 

space of A. The need often arises however for the orthogonal projection onto some arbitrary subspace M. 

The key to using the old P is simply to realize that every subspace is the column space of some matrix. 

More precisely, if 

{xi,...,x m } (5.19) 

is a basis for M then clearly if these Xj are placed into the columns of a matrix called A then K (A) = M. 

For example, if M is the line through ( 1 1 ) then 



=T 



is 




p = \ U ( J. i J 

(5.20) 



is orthogonal projection onto M. 
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5.1.5 Exercises 

1. Gilbert Strang was stretched on a rack to lengths £ = 6, 7, and 8 feet under applied forces of / = 1, 2, 
and 4 tons. Assuming Hooke's law £ — L = cf, find his compliance, c, and original height, L, by least 
squares. 

2. With regard to the example of §3 note that, due to the the random generation of the noise that taints 
the displacements, one gets a different 'answer' every time the code is invoked. 

a. Write a loop that invokes the code a statistically significant number of times and submit bar plots 
of the average fiber stiffness and its standard deviation for each fiber, along with the associated 
M-file. 

b. Experiment with various noise levels with the goal of determining the level above which it becomes 
difficult to discern the stiff fiber. Carefully explain your findings. 

3. Find the matrix that projects K 3 onto the line spanned by I 1 1 ) • 

4. Find the matrix that projects M 3 onto the plane spanned by I 1 1 ) and 1 1 — 1 ) . 

5. If P is the projection of M m onto a fc-dimensional subspace M, what is the rank of P and what is 
R(P)? 



Chapter 6 

Matrix Methods for Dynamical Systems 



6.1 Nerve Fibers and the Dynamic Strang Quartet 1 

6.1.1 Introduction 

Up to this point we have largely been concerned with 

1. Deriving linear systems of algebraic equations (from considerations of static equilibrium) and 

2. The solution of such systems via Gaussian elimination. 

In this module we hope to begin to persuade the reader that our tools extend in a natural fashion to the 
class of dynamic processes. More precisely, we shall argue that 

1. Matrix Algebra plays a central role in the derivation of mathematical models of dynamical systems 
and that, 

2. With the aid of the Laplace transform in an analytical setting or the Backward Euler method in the 
numerical setting, Gaussian elimination indeed produces the solution. 

6.1.2 Nerve Fibers and the Dynamic Strang Quartet 

6.1.2.1 Gathering Information 

A nerve fiber's natural electrical stimulus is not direct current but rather a short burst of current, the so- 
called nervous impulse. In such a dynamic environment the cell's membrane behaves not only like a leaky 
conductor but also like a charge separator, or capacitor. 



1 This content is available online at <http://cnx.Org/content/ml0168/2.6/>. 
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An RC model of a nerve fiber 




Figure 6.1 



The typical value of a cell's membrane capacitance is 

C= 1 2 

cm z 

where /zF denotes micro-Farad. Recalling our variable conventions (Section 2.1.1: Nerve Fibers and the 
Strang Quartet), the capacitance of a single compartment is 

Cm = 2ira — c 

N 

and runs parallel to each R m , see Figure 6.1 (An RC model of a nerve fiber). This figure also differs from 
the simpler circuit (Figure 2.3: The fully dressed circuit model) from the introductory electrical modeling 
module in that it possesses two edges to the left of the stimuli. These edges serve to mimic that portion of 
the stimulus current that is shunted by the cell body. If A c ^ denotes the surface area of the cell body, then 
it has 

Definition 6.1: capacitance of cell body 

L^cb = A c foC 

Definition 6.2: resistance of cell body 

-Kcb = -^cbPrn- 



6.1.2.2 Updating the Strang Quartet 

We ask now how the static Strang Quartet (p. 8) of the introductory electrical module should be augmented. 

6.1.2.2.1 Updating (SI') 

Regarding (SI') (p. 13) we proceed as before. The voltage drops are 

ei = xi 
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e 2 = x\ - E„ 



e 3 = xi - x 2 



e 4 = x 2 



e-h = x 2 - E„ 



e e = x 2 - x 3 



e 7 = x 3 



x 3 - E„ 



and so 



Ax where b = (—E m ) 



( °\ 


1 








1 








V 1/ 



and 



1 -1 








-1 








-1 


1 








-1 








-1 








-1 


1 










V o 








!/ 



6.1.2.2.2 Updating (S2) 

To update (S2) (Section 2.1.2.2: Strang Quartet, Step 2) we must now augment Ohm's law with 

Definition 6.3: Voltage-current law obeyed by a capacitor 

The current through a capacitor is proportional to the time rate of change of the potential across 



it. 



This yields, (denoting derivative by '), 



2/i = C ch ei 



<2 



y^ 


Rcb 


V3 


e 3 

Ri 


2/4 = 


C m e4 


2/5 = 


e 5 




e& 



2/o 



Ri 
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2/7 = C m e 7 ' 



or, in matrix terms. 



where 



and 



2/8 



R n 



Ge + Ce! 
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^m 
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0/ 



are the conductance and capacitance matrices. 



6.1.2.2.3 Updating (S3) 



As Kirchhoff's Current law is insensitive to the type of device occupying an edge, step (S3) proceeds exactly 
as before (Section 2.1.2.3: Strang Quartet, Step 3). 

if) ~ 2/1 - 2/2 - 2/3 = 



2/3 - 2/4 - 2/5 -2/6 = 



or, in matrix terms, 



2/6-2/7-2/8 = 



Ay = — / where / 



( l °\ 
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6.1.2.2.4 Step (S4): Assembling 

Step (S4) remains one of assembling, 

(A T y = -f) =► (A T (Ge + Ce') = -/) =» (A T (G (b - Ax) + C {b> - Ax')) = -/) 

becomes 

A T CAx' + A T GAx = A T Gb + f + A T Cb'. 



(6.1) 



This is the general form of the potential equations for an RC circuit. It presumes of the user knowledge 
of the initial value of each of the potentials, 

x(0) = X (6.2) 

Regarding the circuit of Figure 6.1 (An RC model of a nerve fiber), and letting G = i, we find 

( c cb o o\ 

C 



A T CA 



A T GA 



\ c J 



A T Gb = E n 
and an initial (rest) potential of 



/ G cb + G % -Gi \ 

— Gi 2Gi + G m —Gi 
\ — Gi Gi + G m ) 

( G ch \ { o\ 

and A T Cb' = 



G m 
\ G m j 



x (0) = E„ 





v o y 



i 
v i/ 



6.1.2.3 Modes of Attack 

We shall now outline two modes of attack on such problems. The Laplace Transform (Section 6.2) is an 
analytical tool that produces exact, closed-form solutions for small tractable systems and therefore offers 
insight into how larger systems 'should' behave. The Backward-Euler method (Section 6.4) is a technique 
for solving a discretized (and therefore approximate) version of (6.1). It is highly flexible, easy to code, and 
works on problems of great size. Both the Backward-Euler and Laplace Transform methods require, at their 
core, the algebraic solution of a linear system of equations. In deriving these methods we shall find it more 
convenient to proceed from the generic system 

x' = Bx + g (6.3) 

With respect to our fiber problem 

B = (-(A T CA)~ 1 )A T GA 

\ 

(6.4) 



and 



/ 


-(Gcb+Gi) 


Gi 





\ 




C cb 


C cb 






Gi 


~(2Gi+G m ) 


G, 






c„, 


C m 


Cm 




V 





Gi 

c„, 


-(Gi + G m 
C m 


l ) 



g=(A T CA) 1 (A T Gb+f) 



/ G c b-E m +in \ 



C cb 

E m G„ 

C m 
E m G„ 
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6.2 The Laplace Transform 2 

The Laplace Transform is typically credited with taking dynamical problems into static problems. Recall 
that the Laplace Transform of the function h is 

/>oo 

C{h{s))= / e- {st) h{t)dt. 
Jo 

MATLAB is very adept at such things. For example: 

Example 6.1: The Laplace Transform in MATLAB 

3> syms t 

3> laplace(exp(t) ) 

ans = l/(s-l) 

3> laplace(t* (exp(-t) ) 

ans = l/(s+l)~2 

The Laplace Transform of a matrix of functions is simply the matrix of Laplace transforms of the individual 
elements. 

Example 6.2: Laplace Transform of a matrix of functions 



te 



i 

s-l 



Now, in preparing to apply the Laplace transform to our equation from the dynamic Strang quartet module 
(6.3): 

x' = Bx + g, 

we write it as 

C^j=C(Bx + g) (6.5) 

and so must determine how C acts on derivatives and sums. With respect to the latter it follows directly 
from the definition that 

C{Bx + g) = C {Bx) + C (g) 

= BC(x)+C(g) ' 
Regarding its effect on the derivative we find, on integrating by parts, that 

C (*f) = j°° e-(««)M)(fc = x (t) e-( s ')|S° + a [°° ^ (st) ^ (t) dt. 
\dtj 7o dt ' J 

Supposing that x and s are such that x (t) e~( st ^ — > as t — » oo we arrive at 

c( d ^\=sC(x)-x(Q). (6.7) 



(6.6) 



2 This content is available online at <http://cnx.Org/content/ml0169/2.5/>. 
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Now, upon substituting (6.6) and (6.7) into (6.5) we find 

sC (x) -x(0) = BC (x) + C(g), 
which is easily recognized to be a linear system for C (x), namely 

(si - B)C(x) = C (g) + x (0) . 



(6i 



The only thing that distinguishes this system from those encountered since our first brush (Section 2.1) with 
these systems is the presence of the complex variable s. This complicates the mechanical steps of Gaussian 
Elimination or the Gauss- Jordan Method but the methods indeed apply without change. Taking up the 
latter method, we write 

jC(x) = (sI-B)' 1 (C(g)+x(0)). 

The matrix (si — B)~~ is typically called the transfer function or resolvent, associated with B, at s. 
We turn to MATLAB for its symbolic calculation, (for more information, see the tutorial 3 on MATLAB's 
symbolic toolbox). For example, 

Example 6.3 

> B = [2 -1; -1 2] 

» R = inv(s*eye(2)-B) 

R = 

[ (s-2)/(s*s-4*s+3), -l/(s*s-4*s+3)] 

[ -l/(s*s-4*s+3) , (s-2)/(s*s-4*s+3)] 



We note that (si — B)~ is well defined except at the roots of the quadratic, s 2 — 4s + 3. This quadratic is 
the determinant of si — B and is often referred to as the characteristic polynomial of B. Its roots are 
called the eigenvalues of B. 

Example 6.4 

As a second example let us take the B matrix of the dynamic Strang quartet module (6.4) with 
the parameter choices specified in fib3.m 4 , namely 

/ -0.135 0.125 
0.5 -1.01 



B 



\ 



\ 







0.5 





0.5 
-0.51 / 



(6.9) 



The associated (si — B) is a bit bulky (please run fib3.m 5 ) so we display here only the denom- 
inator of each term, i.e., 



s 6 + 1.655s 2 + 0.4078s + 0.0039. 

3 -A 

Assuming a current stimulus of the form Iq (t) = % n e » ° and E m = brings 



(6.10) 



C(9)(*) 



tee j tp 

10000 *" m - C/ " 

/ 0.191 \ 

(sHY 

o 



V 



o 



/ 



3 http 
4 http 
5 http 



// www.mathworks.com/access/helpdesk/help/toolbox/symbolic/symbolic.shtml 
// www.caam.rice.edu/-~caam335/cox/lectures/fib3.rn 
// www.caam.rice.edu/~caam335/cox/lectures/fib3.rn 
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and so (6.10) persists in 

( s 2 + 1.5s + 0.27 \ 
C (x) = (si - B)- l L (g) = - — 4 — 



(s + |) (s 3 + 1.655s 2 + 0.4078s + 0.0039) 



0.5s + 0.26 
\ 0.2497 ) 



Now comes the rub. A simple linear solve (or inversion) has left us with the Laplace transform of x. The 
accursed 

Theorem 6.1: No Free Lunch Theorem 

We shall have to do some work in order to recover x from C (x). 

confronts us. We shall face it down in the Inverse Laplace module (Section 6.3). 

6.3 The Inverse Laplace Transform 6 

6.3.1 To Come 

In The Transfer Function (Section 9.2) we shall establish that the inverse Laplace transform of a function h 

is 

1 l'°° 
C- 1 (h)(t) = — e^ + y^h((c+yi)t)dy (6.11) 

^ J-oo 

where i = y/2 — 1 and the real number c is chosen so that all of the singularities of h lie to the left of the 
line of integration. 

6.3.2 Proceeding with the Inverse Laplace Transform 

With the inverse Laplace transform one may express the solution of x' = Bx + g , as 

x (t) = /T 1 ((s/ - By 1 ) (C (g) + x (0)) (6.12) 

As an example, let us take the first component of C (x), namely 

0.19 (s 2 + 1.5s + 0.27) 

An (s) = t • 

(s + |) (s 3 + 1.655s 2 + 0.4078s + 0.0039) 

We define: 

Definition 6.4: poles 

Also called singularities, these are the points s at which C Xl (s) blows up. 

These are clearly the roots of its denominator, namely 

»/273 

-1/100, -329/400+^ , and -1/6. (6.13) 

16 

All four being negative, it suffices to take c = and so the integration in (6.11) proceeds up the imaginary 
axis. We don't suppose the reader to have already encountered integration in the complex plane but hope 
that this example might provide the motivation necessary for a brief overview of such. Before that however 
we note that MATLAB has digested the calculus we wish to develop. Referring again to fib3.m 7 for details 
we note that the ilaplace command produces 



6 This content is available online at <http://cnx.Org/content/ml0170/2.8/>. 
7 http:// www.caam.rice.edu/~caam335/cox/lectures/fib3.rn 
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Xl (t) = 211.35eioo 

e^f 1 (262.842cosh (^fp) 



- (0.0554£ 3 + 4.5464£ 2 + 1.085* + 474.19) e 
262.836sinh (^fp) 




Figure 6.2: The 3 potentials associated with the RC circuit model figure (Figure 6.1: An RC model 
of a nerve fiber). 



The other potentials, see the figure above, possess similar expressions. Please note that each of the poles 
of £(xi) appear as exponents in x\ and that the coefficients of the exponentials are polynomials whose 
degrees is determined by the order of the respective pole. 

6.4 The Backward-Euler Method 8 



Where in the Inverse Laplace Transform (Section 6.3) module we tackled the derivative in 

x' = Bx + g, (6.14) 

via an integral transform we pursue in this section a much simpler strategy, namely, replace the derivative 
with a finite difference quotient. That is, one chooses a small dt and 'replaces' (6.14) with 

x(t)-x{t-dt) 



dt 



Bx (i) + g(t). 



(6.15) 



The utility of (6.15) is that it gives a means of solving for x at the present time, t, from the knowledge of 
x in the immediate past, t — dt. For example, as x (0) = x (0) is supposed known we write (6.15) as 



d~t 



B)x(dt) = ^ + g(dt). 



3 This content is available online at <http://cnx.Org/content/ml0171/2.6/>. 
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Solving this for x (dt) we return to (6.15) and find 

L-B^x{2dt) = X -^+g{2dt). 
and solve for x (2dt). The general step from past to present, 

*U*)=(1-b)~ 1 I *®-™ +9V*)) (6.16) 




dt ) V dt 

is repeated until some desired final time, Tdt, is reached. This equation has been implemented in fib3.m 9 
with dt = 1 and B and g as in the dynamic Strang module (6.4). The resulting x ( run fib3.m 10 yourself!) 
is indistinguishable from the plot we obtained (Figure 6.2) in the Inverse Laplace module. 

Comparing the two representations, this equation (6.12) and (6.16), we see that they both produce the 
solution to the general linear system of ordinary equations, see this eqn (6.3), by simply inverting a shifted 
copy of B. The former representation is hard but exact while the latter is easy but approximate. Of course 
we should expect the approximate solution, x , to approach the exact solution, x, as the time step, dt , 
approaches zero. To see this let us return to (6.16) and assume, for now, that g = . In this case, one can 
reverse the above steps and arrive at the representation 

x (jdt) = ({I - dtBY l \ 'x (0) (6.17) 

Now, for a fixed time t we suppose that dt = -. and ask whether 

x (t) = limit 

j— >oo 

This limit, at least when B is one- by-one, yields the exponential 

x(t) = e Bt x(0) 

clearly the correct solution to this equation (6.3). A careful explication of the matrix exponential and its 
relationship to this equation (6.12) will have to wait until we have mastered the inverse laplace transform. 

6.5 Exercises: Matrix Methods for Dynamical Systems 11 

1. Compute, without the aid of a machine, the Laplace transforms of e* and te~ l . Show ALL of your 
work. 

2. Extract from fib3.m analytical expressions for X2 and x%. 

3. Use eig to compute the eigenvalues of B as given in this equation (6.9). Use det to compute the 
characteristic polynomial of B. Use roots to compute the roots of this characteristic polynomial. 
Compare these to the results of eig. How does Matlab compute the roots of a polynomial? (type help 
roots for the answer). 

4. Adapt the Backward Euler portion of fib3.m so that one may specify an arbitrary number of com- 
partments, as in f ibl .m. Submit your well documented M-file along with a plot of x\ and x\$ versus 
time (on the same well labeled graph) for a nine compartment fiber of length I = 1cm. 



9 http:// www.caam.rice.edu/~caam335/cox/lectures/fib3.rn 
10 http:// www.caam.rice.edu/~caam335/cox/lectures/fib3.rn 
11 This content is available online at <http://cnx.Org/content/ml0526/2.4/>. 
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5. Derive this equation (6.17) from a previous equation (6.16) by working backwards toward x (0). Along 
the way you should explain why 



d(t) 



(I-d(t)BY 



6. Show, for scalar B, that \[l-\B 



-i\ J 



e as j — > oo. Hint: By definition 



1- -B 



-iV 



J log 



now use L'Hopital's rule to show that jlog .^ „ — > Bt. 



6.6 Supplemental 

6.6.1 Matrix Analysis of the Branched Dendrite Nerve Fiber 12 

6.6.1.1 Introduction 

In the prior modules on static (Section 2.1) and dynamic electrical systems, we analyzed basic, hypothetical 
one-branch nerve fibers using a modeling methodology we dubbed the Strang Quartet. You may be asking 
yourself whether this method is stout enough to handle the real fiber of our minds. Indeed, can we use our 
tools in a real-world setting? 

6.6.1.2 Presentation 



An Actual Nerve Fiber 




Figure 6.3: A pyramidal neuron from the CA3 region of a rat's hippocampus, scanned at (FIX ME) X 
magnification. 



12 This content is available online at <http://cnx.Org/content/ml0177/2.7/>. 
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To answer your question, the above is a rendering of a neuron from a rat's hippocampus. The tools we 
have refined will enable us to model the electrical properties of a dendrite leaving the neuron's cell body. A 
three-branch model of such a dendrite, traced out with painstaking accuracy, appears in the diagram below. 



3-branch Dendrite Model 







w^w 



WF^WP 



w 



Yij. 



'Jl^i _L ^«_ 




Figure 6.4: Multi-compartment electrical model of a rendered dendrite fiber. 



Our multi-compartment model reveals a 3 branch, 10 node, 27 edge structure to the fiber. Note that 
we have included the Nernst potentials (Definition: "Nernst potentials", p. 12), the nervous impulse as a 
current source, and the additional leftmost edges depicting stimulus current shunted by the cell body. 

We will continue using our previous notation, namely: Ri and R m denoting cell body (Definition: "axial 
resistance", p. 6) and membrane (Definition: "membrane resistance", p. 6) resistances, respectively; x 



representing the vector of potentials x\ . . . Xio, and y denoting the vector of currents y\ . . 
typical value for a cell's membrane capacitance, 

c= l(fiF/cm 2 ), 

we derive (see variable conventions (Section 2.1.1: Nerve Fibers and the Strang Quartet)): 
Definition 6.5: Capacitance of a Single Compartment 



. 7/27- Using the 



G„ 



I 

2ira—c 

N 



This capacitance is modeled in parallel with the cell's membrane resistance. Additionally, letting A c \~, 
denote the cell body's surface area, we recall that its capacitance and resistance are 

Definition 6.6: Capacitance of cell body 

L'cb = A C \~,C 



65 
Definition 6.7: Resistance of cell body 

6.6.1.3 Applying the Strang Quartet 

6.6.1.3.1 Step (Sl')-Voltage Drops 

Let's begin filling out the Strang Quartet. For Step (SI'), we first observe the voltage drops in the figure. 
Since there are a whopping 27 of them, we include only the first six, which are slightly more than we need 
to cover all variations in the set: 

ei = xi 



e 2 = xi- E n 



e 3 = xi - x 2 



e 4 = x 2 



es = x 2 - E n 



e 6 = x 2 - x 3 . 



e 2 7 = X\o — E„ 
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In matrix for, letting b denote the vector of batteries, 



b — Ax where 



{—E m ) 



o\ 

1 




1 




1 





1 




1 




1 




1 




1 
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and 
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Although our adjacency matrix A is appreciably larger than our previous examples, we have captured the 
same phenomena as before. 



6.6.1.3.2 Applying (S2): Ohm's Law Augmented with Volt age- Current Law for Capacitors 

Now, recalling Ohm's Law and remembering that the current through a capacitor varies proportionately 
with the time rate of change of the potential across it, we assemble our vector of currents. As before, we list 
only enough of the 27 currents to fully characterize the set: 



I/i= C, 



cb 



~di 



1)2 



Rcb 
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In matrix terms, this compiles to 



where 



2/3 


e 3 

Ri 


2/4 = 


de 4 


2/5 = 


e 5 


2/27 


e 27 

■Km 


y=G 





Conductance matrix 
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(6.18) 
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Capacitance matrix 
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6.6.1.3.3 Step (S3): Applying Kirchoff's Law 

Our next step is to write out the equations for Kirchoff's Current Law. We see: 

io ~ 2/1 - 2/2 - 2/3 = 

2/3-2/4-2/5-2/6 = 

2/6-2/7-2/8-2/9 = 

2/9 - 2/io - 2/19 = 

2/io - 2/ii - 2/12 - 2/13 = 

2/io - 2/ii - 2/12 - 2/13 = 

2/13 - 2/14 - 2/15 - 2/16 = 

2/16 - 2/17 - 2/18 = 

2/19 - 2/20 - 2/21 - 2/22 = 

2/22 - 2/23 - 2/24 - 2/25 = 

2/25 - 2/26 - 2/27 = 

Since the B coefficient matrix we'd form here is equal to A T , we can say in matrix terms: 

A T y = -f 
where the vector / is composed of /i = Iq and /2...27 = 0. 

6.6.1.3.4 Step (S4): Stirring the Ingredients Together 

Step (S4) directs us to assemble our previous toils together into a final equation, which we will then endeavor 
to solve. Using the process (Section 6.1.2.2.4: Step (S4): Assembling) derived in the dynamic Strang module, 
we arrive at the equation 

A T CA— + A T GAx = A T Gb+f + A T C-, (6.20) 

at at 

which is the general form for RC circuit potential equations. As we have mentioned, this equation presumes 

knowledge of the initial value of each of the potentials, x (0) = X. 



72 



CHAPTER 6. MATRIX METHODS FOR DYNAMICAL SYSTEMS 



Observing our circuit (Figure 6.4: 3-branch Dendrite Model), and letting -^-j- 
necessary quantities to fill out (6.20) 's pieces (for these calculations, see dendrite. m 



Gfoo! we calculate the 
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5 http://www.ece. rice.edu/-~rainking/dendrite/matlab/dendrite.rn 
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and an initial (rest) potential of 



A 1 C 



db 
~dt 



x (0) = E„ 



0, 

( l\ 

1 
1 
1 
1 
1 
1 
1 
1 

v 1 ; 



6.6.1.4 Applying the Backward-Euler Method 

Since our system is so large, the Backward-Euler method is the best path to a solution. Looking at the 
matrix A T CA, we observe that it is singular and therefore non-invertible. This singularity arises from the 
node connecting the three branches of the fiber and prevents us from using the simple equation x' = Bx + g, 
we used in earlier Backward-Euler-ings (Section 6.4). However, we will see that a modest generalization to 
our previous form yields (6.21): 

Dx' = Ex + g (6.21) 

capturing the form of our system and allowing us to solve for x (t). We manipulate (6.21) as follows: 

Dx' = Ex+ q 



x(i)-x(t- dt) 



where in our case 



(D - Edt) x (i) = Dx(t- dt) + gdt 



x (t) = (D - Edt) (Dx (t - dt) + gdt) , 



D = A T CA, 



E = - (A T GA) , and 



g = A T Gb+f. 

This method is implemented in dendrite. m 14 with typical cell dimensions and resistivity properties, yielding 
the following graph of potentials. 



* http://it.is. rice.edu/-~rainking/dendrite/matlab/dendrite.iri 
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Graph of Dendrite Potentials 



E 




(a) 




M 35 



Time (msec] 



(b) 



Figure 6.5: (a) Large view of potentials, (b) Zoomed view of potentials showing maxima. 



Chapter 7 

Complex Analysis 1 



7.1 Complex Numbers, Vectors and Matrices 1 

7.1.1 Complex Numbers 

A complex number is simply a pair of real numbers. In order to stress however that the two arithmetics 
differ we separate the two real pieces by the symbol i. More precisely, each complex number, z, may be 
uniquely expressed by the combination x + iy, where x and y are real and i denotes y/—l. We call x the real 
part and y the imaginary part of z. We now summarize the main rules of complex arithmetic. 
If Z\ = X\ + iy\ and z 2 = x 2 + iy 2 then 
Definition 7.1: Complex Addition 
z\ + z 2 = xi + x 2 + i (yi + 2/2) 
Definition 7.2: Complex Multiplication 
ziz 2 = (xi + iyi) (x 2 + iy 2 ) = xix 2 - yiy 2 + i (x\y 2 + x 2 y\) 
Definition 7.3: Complex Conjugation 
zl = x\- iyx 
Definition 7.4: Complex Division 

£1 _ £1 zj _ x 1 x 2 +yiy2+i(x2Vi-xiy 2 ) 
Z 2 ~ Z 2 Z2 — x 2 ' 2 +y 2 2 

Definition 7.5: Magnitude of a Complex Number 

\zi\ = \Jz{zi= -Jx-y 2 + yx 2 

7.1.2 Polar Representation 

In addition to the Cartesian representation z = x + iy one also has the polar form 

z = \z\(cos{6) + ism{6)) (7.1) 

where 8 = arctan (-). 

This form is especially convenient with regards to multiplication. More precisely, 

ziz 2 = \z 1 \\z 2 \(cos(9i)cos(9 2 ) -s'm(8 1 )sm(6 2 ) + i(cos(8 1 )sm(9 2 ) + sm(6 1 )cos(6 2 ))) 
= |2i||2 2 |(cos(0i+02)+«sin(0 1 + 2 )) 

As a result: 

z n = (\ z \) n (cos (n6) + isin (n9)) (7.3) 



1 This content is available online at <http://cnx.Org/content/ml0504/2.5/>. 
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7.1.3 Complex Vectors and Matrices 

A complex vector (matrix) is simply a vector (matrix) of complex numbers. Vector and matrix addition 
proceed, as in the real case, from elementwise addition. The dot or inner product of two complex vectors 
requires, however, a little modification. This is evident when we try to use the old notion to define the length 
of a complex vector. To wit, note that if: 

I 1 + i \ 



then 



T 

z z 



(i + iY + (i-iY 



i 



l + 2z— 1 + 1 — 2f— 1 = 



Now length should measure the distance from a point to the origin and should only be zero for the zero 
vector. The fix, as you have probably guessed, is to sum the squares of the magnitudes of the components 
of z. This is accomplished by simply conjugating one of the vectors. Namely, we define the length of a 
complex vector via: 

(z) = "JWz (7.4) 



In the example above this produces 



ll + il 



\l-i\f = y/l=2 



As each real number is the conjugate of itself, this new definition subsumes its real counterpart. 

The notion of magnitude also gives us a way to define limits and hence will permit us to introduce 

( M 

complex calculus. We say that the sequence of complex numbers, < 



z n \n 



complex number zq and write 



2 

V ■■■ ) 



>, converges to the 



Z n —* Zq 



Zq = limit z n 

n — *oo 

when, presented with any e > one can produce an integer N for which \z n — zq\ < e when n > N. As an 
example, we note that (|) — * 0. 

7.1.4 Examples 

Example 7.1 

As an example both of a complex matrix and some of the rules of complex arithmetic, let us 
examine the following matrix: 



(l 


1 


1 


1 \ 


1 


i 


-1 


— i 


1 

1 


-1 


1 
1 


-1 



V ! 



1 



(7.5) 



/ 



Let us attempt to find FF. One option is simply to multiply the two matrices by brute force, 
but this particular matrix has some remarkable qualities that make the job significantly easier. 
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Specifically, we can note that every element not on the diagonal of the resultant matrix is equal to 
0. Furthermore, each element on the diagonal is 4. Hence, we quickly arrive at the matrix 

\ 



FF 



1 4 














4 














4 





V° 








4 



(7.6) 



= U 

This final observation, that this matrix multiplied by its transpose yields a constant times the 
identity matrix, is indeed remarkable. This particular matrix is an example of a Fourier matrix, 
and enjoys a number of interesting properties. The property outlined above can be generalized for 
any F n , where F refers to a Fourier matrix with n rows and columns: 



F n F n = nl 



(7.7) 



7.2 Complex Functions 2 
7.2.1 Complex Functions 

A complex function is merely a rule for assigning certain complex numbers to other complex numbers. The 
simplest (nonconstant) assignment is the identity function / (z) = z. Perhaps the next simplest function 
assigns to each number its square, i.e., / (z) = z 2 . As we decomposed the argument of /, namely z, into 
its real and imaginary parts, we shall also find it convenient to partition the value of /, z 2 in this case, into 
its real and imaginary parts. In general, we write 

f(x + iy) = u(x,y)+iv(x,y) (7.8) 

where u and v are both real- valued functions of two real variables. In the case that / (z) = z 2 we find 

u(x,y) =x 2 - y 2 

and 

v(x,y) = 2xy 

With the tools of complex numbers (Section 7.1), we may produce complex polynomials 

/ (z) = z m + c m _ 1 2 m - 1 + • • • + c t z + c (7.9) 

We say that such an / is order m. We shall often find it convenient to represent polynomials as the product 
of their factors, namely 



/(z) = (z-A 1 ) dl (z-A 2 ) d2 



{z-\ h ) ah 



(7.10) 



Each Xj is a root of degree dj. Here h is the number of distinct roots of /. We call Xj a simple root when 
dj = 1. We can observe the appearance of ratios of polynomials or so called rational functions. Suppose 



«(*) 



/(*) 



2 This content is available online at <http://cnx.Org/content/ml0505/2.7/>. 
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in rational, that / is of order at most m — 1 while g is of order m with the simple roots {Ai, . . . , A m }. It 
should come as no surprise that such q should admit a Partial Fraction Expansion 

m 

One uncovers the g.,- by first multiplying each side by z — Xj and then letting z tend to Xj. For example, if 

z 2 + 1 2: + i z — i 
then multiplying each side by z + i produces 



(7.12) 



— = ft + ^±^ (7.13) 



Now, in order to isolate (ft it is clear that we should set z = —i. So doing we find that q^ = |. In order to 
find q 2 we multiply (7.12) by z — i. and then set z = i. So doing we find q 2 = ^, and so 



1 



z 2 + i z + i z — i 
. Returning to the general case, we encode the above in the simple formula 



(7.14) 



qj = limit (z — Xj) q (z) (7-15) 

zz — >A, 



You should be able to use this to confirm that 

z 1/2 1/2 



z 2 + 1 z + i 



(7.16) 



Recall that the transfer function we met in The Laplace Transform (p. 58) module was in fact a matrix 
of rational functions. Now, the partial fraction expansion of a matrix of rational functions is simply the 
matrix of partial fraction expansions of each of its elements. This is easier done than said. For example, the 
transfer function of 

. 1 
B 

-1 



(zI-B)- 1 




(7.17) 



The first line comes form either Gauss- Jordan by hand or via the symbolic toolbox in Matlab. More 
importantly, the second line is simply an amalgamation of (7.12) and (7.14). Complex matrices have finally 
entered the picture. We shall devote all of Chapter 10 to uncovering the remarkable properties enjoyed 
by the matrices that appear in the partial fraction expansion of (zl — B)~ Have you noticed that, in our 
example, the two matrices are each projections, and they sum to /, and that their product is 0? Could this 
be an accident? 

In The Laplace Transform (Section 6.2) module we were confronted with the complex exponential. By 
analogy to the real exponential we define 

oo n 

e* = £ Z - (7.18) 

n=0 
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and find that 

e e = 1 + ^+M)! + M! + M! + ... 

= i-^ + S---. + <(tf-ft + ft-...) (7.19) 

= cos (0) + isin (8) 

With this observation, the polar form (Section 7.1.2: Polar Representation) is now simply z = \z\e t0 . 

One may just as easily verify that 

i9 , J-i)e 

cos{9) = —^ 

and 

id _ J-i)9 

sin (0) = 



% 

These suggest the definitions, for complex z, of 



e iz _|_ M-z) 

cos(z) = (7.20) 



and 

e iz — e (-i) z 



sin(z) = (7.21) 

2t 



As in the real case the exponential enjoys the property that 



e zi+z 2 = e z 1& z 2 



and in particular 

p x+iy 



f(Zn) ~ f(zo) 

same value for every seque 

the limit £f{z ) 

3 This content is available online at <http://cnx.Org/content/ml0276/2.7/>. 



(7.22) 



= e^cos (y) + ie^sin (y) 
Finally, the inverse of the complex exponential is the complex logarithm, 

]n(z)=hi(\z\) + i9 (7.23) 

for z = \z\e 10 . One finds that In (-1 + i) = In (a/2) + i^f. 

7.3 Complex Differentiation 3 

7.3.1 Complex Differentiation 

The complex / is said to be differentiable at zq if 

limit '('WW 

z— >z Z — Zq 

exists, by which we mean that 



z n — Zq 
converges to the same value for every sequence {z n } that converges to zq. In this case we naturally call 
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To illustrate the concept of 'for every' mentioned above, we utilize the following picture. We assume 
the point zq is differentiable, which means that any conceivable sequence is going to converge to zq. We 
outline three sequences in the picture: real numbers, imaginary numbers, and a spiral pattern of both. 



Sequences Approaching A Point In The Complex Plane 




Figure 7.1: The green is real, the blue is imaginary, and the red is the spiral. 



81 



7.3.2 Examples 



Example 7.2 






The derivative of z 2 is 2z. 








limit 7 - 2 ~ z ° 2 -- 


= limit (*-^+^ 

z^>z z ~ z ° 



2z 



(7.24) 



Example 7.3 

The exponential is its own derivative. 



e*-*0-l 



limit - — - — = e z ° limit 
z—>z z ~ z ° z^z a z ~ z ° 



limit E^o^iyr ( 7 - 25 ) 



e 

z— >z 



= e z ° 

Example 7.4 

The real part of z is not a differentiable function of z. 

We show that the limit depends on the angle of approach. First, when z n — > z$ on a line parallel 
to the real axis, e.g., z n = xq -\ \- iyo, we find 

limit *°+*-*° = ! (7 . 26) 

™^°o x + - + ij/o - x + iyo 



while if z n — > zq in the imaginary direction, e.g., z n = xq + i (yo + ), then 



limit , X ° .f = (7.27) 

- cc + i (2/0 + ;J - x + iy 



n —> oc 



7.3.3 Conclusion 

This last example suggests that when / is differentiable a simple relationship must bind its partial derivatives 
in x and y. 

Proposition 7.1: Partial Derivative Relationship 

If / is differentiable at z then fj (z ) = ^^ = - (i^ 1 ) 



Proof: 

With z = x + iyo, 



£/(*>) = " ll ^ A 



Z—>Zo 

limit /( a: +^)-/( a; °+^o) (7.28) 

x— >xo x x ° 

df(zp) 

dx 



Alternatively, when z = Xq 4- iy then 

4-f(zo) = limit /(z) - /(zo) 



z-zo 



z— >z 

/(xo+ii/)-/(xQ+ii/o) 

y^yo 



limit /('°+'?/)-/(^+'«°) (7.29) 



O 2 ^) 
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7.3.4 Cauchy-Reimann Equations 

In terms of the real and imaginary parts of / this result brings the Cauchy-Riemann equations. 

du dv 
dx dy 

and 

Regarding the converse proposition we note that when / has continuous partial derivatives in region obeying 
the Cauchy-Reimann equations then / is in fact differentiable in the region. 

We remark that with no more energy than that expended on their real cousins one may uncover the rules 
for differentiating complex sums, products, quotients, and compositions. 

As one important application of the derivative let us attempt to expand in partial fractions a rational 
function whose denominator has a root with degree larger than one. As a warm-up let us try to find gii 
and q\ t 2 in the expression 

z+2 _ (? M ffl,2 

{z+lf ~ z+1 {z+lf 

Arguing as above, it seems wise to multiply through by (z + 1) and so arrive at 

z + 2 = q 1A (z + l) + q h 2 (7.32) 

On setting z = —\ this gives q 12 = 1. With q 12 computed, (7.32) takes the simple form z + 1 = q\\ (z + 1) 
and so gi 2 = 1 as well. Hence, 

z+2 _ 1 1 

{z + lf ~ ~z~+l{z+lf 

This latter step grows more cumbersome for roots of higher degrees. Let us consider 

{z + 2) 2 _ q lfl q 1>2 qi,3 



{z+lf z+1 {z+lf {z+lf 

The first step is still correct: multiply through by the factor at its highest degree, here 3. This leaves us 
with 

{z + 2) 2 = qi>1 (z + if + qi,2 (Z + 1)+ qi, 3 (7.33) 

Setting z = —1 again produces the last coefficient, here qi 3 = 1. We are left however with one equation in 
two unknowns. Well, not really one equation, for (7.33) is to hold for all z. We exploit this by taking two 
derivatives, with respect to z, of (7.33). This produces 

2{z + 2) = 2q 1 . 1 {z+l)+q h2 

and 2 = gr^i The latter of course needs no comment. We derive q\ t2 from the former by setting z = —1. This 
example will permit us to derive a simple expression for the partial fraction expansion of the general 
proper rational function q = *- where g has h distinct roots {Ai, . . . , A^} of respective degrees {di, . . . , d^}. 
We write 

ffC^EEr^Tv* (7 - 34) 
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d:—k . 



and note, as above, that qj^ is the coefficient of (z — dj) 3 in the rational function 

rj(z) = q(z)(z-\ J ) dl 

Hence, q j j, may be computed by setting z = Xj in the ratio of the dj — fcth derivative of tj to (dj — k)l. 
That is, 

q jik = limit — ^ ,,^,-k \( z ~ X jf j< l( z )} ( 7 - 35 ) 



z^\j (dj - k)\ dz° 



As a second example, let us take 



B 



( 1 \ 
1 3 

\ o i i y 



(7.36) 



and compute the $,• h matrices in the expansion 



(zI-B)- 1 



( ^ 



\ 



\ 



(z-l) 2 (z-3) (z-l)(z-3) 7^1 I 



(z-l)(z-3) 
1 




1 

z-3 

1 



Z- 1 



*1, 



(z-iy 



r$l,: 



z-3 



$2,1 



The only challenging term is the (3, 1) element. We write 

1 9i,i 9i,2 



92,1 



(z-lY(z-S) 



It follows from (7.35) that 



and 



9i,i 



-i (z-iy z-z 

dz \ z— 3 / 

-1/4 



91,2 



and 



It now follows that 



92,1 



^3-1 

-1/4 



z-3 

1/4 



(7.37) 



(7.38) 



(7.39) 



(zI-B)- 1 



z- 1 



' 1 0^ 

-1/2 

V -1/4 "1/2 1 / 



(*-ir 



/ \ 


\ -1/2 0/ 



z-3 



/ \ 

1/2 1 
\ 1/4 1/2 / 



(7.40) 



In closing, let us remark that the method of partial fraction expansions has been implemented in Matlab. 
In fact, (7.37), (7.38), and (7.39) all follow from the single command: [r,p,k]=residue( [0 1] , [1 -5 
7 -3]). The first input argument is Matlab-speak for the polynomial / (z) = 1 while the second argument 
corresponds to the denominator 



g (z) = (z- if (z - 3) = z 3 - bz 2 + 7z - 3 
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7.4 Exercises: Complex Numbers, Vectors, and Functions 4 
7.4,1 Exercises 

Exercise 7.1 (Solution on p. 85.) 

Express \e z \ in terms of x and/or y. 

Exercise 7.2 (Solution on p. 85.) 

Confirm that e ln ^ = z and In (e z ) = z. 

Exercise 7.3 (Solution on p. 85.) 

Find the real and imaginary parts of cos (z) and sin (z). Express your answers in terms of regular 
and hyperbolic trigonometric functions. 

Exercise 7.4 (Solution on p. 85.) 

Show that cos 2 (z) + sin 2 (z) = 1 

Exercise 7.5 (Solution on p. 85.) 

With z w = e win ^ z ' for complex z and w compute \/i 

Exercise 7.6 (Solution on p. 85.) 

Verify that cos (z) and sin (z) satisfy the Cauchy-Riemann equations and use the proposition to 
evaluate their derivatives. 

Exercise 7.7 (Solution on p. 85.) 

Submit a Matlab diary documenting your use of residue in the partial fraction expansion of the 
transfer function of 



B 



( 2 \ 

-14 
V -12/ 



4 This content is available online at <http://cnx.Org/content/ml0506/2.5/>. 
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Solutions to Exercises in Chapter 7 

Solution to Exercise 7.1 (p. 84) 

Pending completion of assignment. 
Solution to Exercise 7.2 (p. 84) 

Pending completion of assignment. 
Solution to Exercise 7.3 (p. 84) 

Pending completion of assignment. 
Solution to Exercise 7.4 (p. 84) 

Pending completion of assignment. 
Solution to Exercise 7.5 (p. 84) 

Pending completion of assignment. 
Solution to Exercise 7.6 (p. 84) 

Pending completion of assignment. 
Solution to Exercise 7.7 (p. 84) 

Pending completion of assignment. 
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Chapter 8 

Complex Analysis 2 



8.1 Cauchy's Theorem 1 

8.1.1 Introduction 

Our main goal is a better understanding of the partial fraction expansion of a given transfer function. With 
respect to the example that closed the discussion of complex differentiation, see this equation (7.36) - In this 
equation (7.40), we found 

{zI ~ Byl = v\ Pl + T^V^ Dl + v\ P2 

z — M (z — Ai) z — a 2 

where the Pj and Dj enjoy the amazing properties 

(8.1) 
= XiPi + D 1 

and 

BP 2 = P 2 B = X 2 P 2 

2. 

P 1 + P 2 = I (8.2) 

Pi 2 = Pi 
P2 2 = P2 



and 
3. 



L>i 2 = 



P 1 D 1 = D 1 P 1 

= D 1 

and 

P 2 D l = D l P 2 = 



*.3) 



In order to show that this always happens, i.e., that it is not a quirk produced by the particular B in this 
equation (7.36), we require a few additional tools from the theory of complex variables. In particular, we 
need the fact that partial fraction expansions may be carried out through complex integration. 



1 This content is available online at <http://cnx.Org/content/ml0264/2.8/>. 
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8.1.2 Integration of Complex Functions Over Complex Curves 

We shall be integrating complex functions over complex curves. Such a curve is parameterized by one 
complex valued or, equivalently, two real valued, function (s) of a real parameter (typically denoted by t). 
More precisely, 

C = {z (t) = x (t) + iy (t) \a<t<b} 

For example, if x (t) = y (t) = t while a = and 6=1, then C is the line segment joining + iO to 1 + i. 
We now define 

f(z)dz = I f(z(t))z'(t)dt 



For example, ifC={£ + i£|0<£<l}as above and f (z) = z then 



zdz= I (t + it)(l + i)dt= / t-t + i2tdt = i 
o Jo 



while if C is the unit circle { e lt | < t < 2ir } then 



27T /.27T />27T 

it-Ati± • I S2t r 



zdz = I e n ie lt dt = i / e lZt dt = i cos (2i) + ism (2i) dt = 

"'O -'o 



Remaining with the unit circle but now integrating f (z) = - we find 

z~ l dz = / e~^\e lt dt = 2m 



We generalize this calculation to arbitrary (integer) powers over arbitrary circles. More precisely, for 
integer m and fixed complex a we integrate (z — a) over 

C{a,r)= {a+ re lt \ < t < 2ir } 

the circle of radius r centered at a. We find 

J (z — a) m dz = L w (a + re %t — a) rie lt dt 



2tt 

I n. -A- rp"" — n. I rrp.'TIt, 

(8.4) 



ir m+1 /g e l(m+1 ^dt 



(z- a) m dz = ir m+1 / cos ((m + l)t) + isin ((m + 1) t) dt 



,; " ! 2m if m = -1 

/o I otherwise 



When integrating more general functions it is often convenient to express the integral in terms of its real 
and imaginary parts. More precisely 



/ (z) dz = / u (x, y) + iv (x, y) dx + i / u (x, y) + iv (x, y) dy 
f (z) dz = / u(x,y)dx— / v(x,y)dy + i / v(x,y)dx + i / u(x,y)dy 



ff(z)dz = f'u(x(t),y(t))x'(t) - v(x(t),y(t))y'(t)dt + 

i j b a u (x (t) , y (*)) y' (t) +v(x(t),y (t)) x f (t) dt 



The second line should invoke memories of: 
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Theorem 8.1: Green's Theorem 

If C is a closed curve and M and N are continuously differentiable real-valued functions on C mj 
the region enclosed by C, then 

f f f dN dM 

Mdx + / Ndy =11- T—dxdy 

J J J ox dy 

Applying this to the situation above, we find, so long as C is closed, that 

. . . , f f dv du , , , f f du dv 

f(z)dz = - / — + T^-dxdy + i / — + -^-dxdy 
J J ox dy J J ox dy 

At first glance it appears that Green's Theorem only serves to muddy the waters. Recalling the Cauchy- 
Riemann equations (Section 7.3.4: Cauchy-Reimann Equations) however we find that each of these double 
integrals is in fact identically zero! In brief, we have proven: 

Theorem 8.2: Cauchy's Theorem 

If / is differentiable on and in the closed curve C then J f (z) dz = 0. 

Strictly speaking, in order to invoke Green's Theorem we require not only that / be differentiable but 
that its derivative in fact be continuous. This however is simply a limitation of our simple mode of proof; 
Cauchy's Theorem is true as stated. 

This theorem, together with (8.4), permits us to integrate every proper rational function. More precisely, 
if q = - where / is a polynomial of degree at most m — 1 and g is an mth degree polynomial with h distinct 
zeros at { Xj \ j = {1, . . . , h} } with respective multiplicities of { rrij \ j = {1, . . . , h} } we found that 

h rrij 

«w = EErfvf (8 - 5) 

Observe now that if we choose rj so small that Xj is the only zero of g encircled by Cj = C (Xj, rj) then by 
Cauchy's Theorem 

q(z)dz = J2lj,k / " —rdz 

k=i J ( Z - X J) 

In (8.4) we found that each, save the first, of the integrals under the sum is in fact zero. Hence, 

q (z) dz = 2niqj t i (8-6) 

With qj t i in hand, say from this equation (7.35) or residue, one may view (8.6) as a means for computing 
the indicated integral. The opposite reading, i.e., that the integral is a convenient means of expressing qj t \, 
will prove just as useful. With that in mind, we note that the remaining residues may be computed as 
integrals of the product of q and the appropriate factor. More precisely, 

q (z) (z — Xj) dz = 2wiqj t k (8-7) 

One may be led to believe that the precision of this result is due to the very special choice of curve and 
function. We shall see ... 

8.2 Cauchy's Integral Formula 2 
8.2.1 The Residue Theorem 

After Cauchy's Theorem eqn6 (8.6) and Cauchy's Theorem eqn7 (8.7) perhaps the most useful consequence 
of Cauchy's Theorem (Theorem 8.2, Cauchy's Theorem, p. 89) is the 



2 This content is available online at <http://cnx.Org/content/ml0246/2.8/>. 
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Lemma 8.1: The Curve Replacement Lemma 

Suppose that C2 is a closed curve that lies inside the region encircled by the closed curve C\. If f 
is differentiable in the annular region outside C2 and inside C\ then 



f(z)dz= / f{z)dz 



Proof: 

With reference to Figure 8.1 (Curve Replacement Figure) we introduce two vertical segments and 
define the closed curves C3 = abcda (where the be arc is clockwise and the da arc is counter- 
clockwise) and C4 = adeba (where the ad arc is counter-clockwise and the cb arc is clockwise). By 
merely following the arrows we learn that 



f(z)dz = J f(z)dz + J f(z)dz + J f(z)dz 



As Cauchy's Theorem (Theorem 8.2, Cauchy's Theorem, p. 89) implies that the integrals over C3 
and C4 each vanish, we have our result. 



Curve Replacement Figure 




Figure 8.1: The Curve Replacement Lemma 
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This Lemma says that in order to integrate a function it suffices to integrate it over regions where it is 
singular, i.e. nondifferentiable. 

Let us apply this reasoning to the integral 

-dz 



(z-\i)(z-X 2 ) 

where C encircles both Ai and A2 as depicted in Figure 8.2. We find that 



z I z I z 

dz = I ; — — ; — ~dz + / ; — — ; — - dz 



(z-X 1 )(z-X 2 ) " J (z-X 1 )(z-X 2 ) " J (z-X 1 )(z-X 2 ) 
Developing the integrand in partial fractions we find 

Ai f 1 A 2 f 1 2niX! 

-dz = — / —dz + — / —dz 



(z — Ai) (z — X 2 ) A x - A 2 J z - Ai A 2 - Ai J z - A 2 Ai - A 2 

Similarly, 

f z _ 27riA 2 

J (z-A!)(z-A 2 ) Z ~ \2^M 
Putting things back together we find 



/ ( z -\ 1 )( z -\ 2 ) dz ~ 2m {\ 1 -X 2 + X 2 -X 1 ) 

= 2iri 



Image notjtnished 

Figure 8.2: Concentrating on the poles. 



We may view (8.8) as a special instance of integrating a rational function around a curve that encircles 
all of the zeros of its denominator. In particular, recalling that cauchy's theorem eqn5 (8.5) and Cauchy's 
theorem eqn6 (8.6), we find 

To take a slightly more complicated example let us integrate ^-^- over some closed curve C inside of 
which / is different iable and a resides. Our Curve Replacement Lemma now permits us to claim that 

m dz= rm dz 



It appears that one can go no further without specifying /. The alert reader however recognizes that in the 
integral over C (a, r) is independent of r and so proceeds to let r — > 0, in which case z — > a and / (z) — > / (a). 
Computing the integral of -^- along the way we are led to the hope that 

^-dz = f (a) 2iri 
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In support of this conclusion we note that 

fMte= [ m + m-m dz=m r j_ dz+ r fw-fw dz 



z — a J z — a z — a z — a J z — a 

Now the first term is / (a) 2iri regardless of r while, as r — > 0, the integrand of the second term approaches 
4-f (a) and the region of integration approaches the point a. Regarding this second term, as the integrand 
remains bounded as the perimeter of C (a, r) approaches zero the value of the integral must itself be zero. 
This result if typically known as 

Rule 8.1: Cauchy's Integral Formula 

If / is differentiable on and in the closed curve C then 

f{a) = i}-.(—d* (8-10) 

2m J z — a 

for each a lying inside C. 

The consequences of such a formula run far and deep. We shall delve into only one or two. First, we note 
that, as a does not lie on C, the right hand side is a perfectly smooth function of a. Hence, differentiating 
each side, we find 

da J 2iriJ (z-a)" 

for each a lying inside C. Applying this reasoning n times we arrive at a formula for the n-th derivative of 

/ at a, 

d n r , x n\ f f(z) , 

f(a) = / V ' dz (8.12) 

da nJ { ' 2m J (z-a) 1+n K ' 

for each a lying inside C. The upshot is that once / is shown to be differentiable it must in fact be infinitely 
differentiable. As a simple extension let us consider 

±r m dz 

2m J ( Z -X l ){z-X 2 f 

where / is still assumed differentiable on and in C and that C encircles both Ai and A2. By the curve 
replacement lemma this integral is the sum 

1 r '(*) cb+J-/- m Kdz 



:/(«) = S3 / T^^dz (8.11) 



2nil (,-\A(?.- \^ 2 2ni 1 * '• ■ ■ ■ ■ ■ :;: 



( z _ Al )(z-A 2 ) z 2m J ( z - Xl )(z-\ 2 y 

where Xj now lies in only Cj. As ^_a ls we ^ behaved in C\ we may use (8.10) to conclude that 

1 /■ Hz) _ M _ /(Ar) 



2 r\ i \2 



2m J ( z -\i)(z-\ 2 ) (A1-A2) 

Similarly, as j^r\ ls we H behaved in C 2 we may use (8.11) to conclude that 

1 f fM 7dz _d(f(a) 



2m J ( z _ X 1 )(z- X 2 f " da\a-Xi, 

These calculations can be read as a concrete instance of 

Theorem 8.3: The Residue Theorem 

If g is a polynomial with roots { Xj \ j = {1, . . . , h} } of degree { dj \ j = {1, . . . , h} } and C is a 
closed curve encircling each of the Xj and / is differentiable on and in C then 

f (z) h 

— -r^rdz = 27ri > res (A 9 ) 

g(z) fri 3 ' 
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where 



res (Xj) = limit 



1 d^'Hz-X^^ 



z^Xj (dj - 1)! dz^- 1 

is called the residue of - at A 7 . 

a J 

One of the most important instances of this theorem is the formula for the Inverse Laplace Transform. 



8.3 The Inverse Laplace Transform: Complex Integration 3 
8.3.1 The Inverse Laplace Transform 

If q is a rational function with poles { Xj \ j = {1, . . . , h} } then the inverse Laplace transform of q is 

£-i( 9 )(t) = -L f ' q { z )e*dz 

2m J 

where C is a curve that encloses each of the poles of q. As a result 

h 
£- 1 (g)(i) = ^res(A j ) 



.13) 



.14) 



Let us put this lovely formula to the test. We take our examples from discussion of the Laplace Transform 
(Section 6.2) and the inverse Laplace Transform (Section 6.3). Let us first compute the inverse Laplace 
Transform of 

q{Z) = ^ 
According to (8.14) it is simply the residue of q (z) e zt at z = — 1, i.e., 

/ x ,- deZt -t 

res (—1) = limit = te 

z—y-l dz 

This closes the circle on the example begun in the discussion of the Laplace Transform (Section 6.2) and 
continued in exercise one for chapter 6. 
For our next example we recall 



C(x 1 (s)) 



0.19 (s 2 + 1.5s + 0.27) 



(s + 1/6) (s 3 + 1.655s 2 + 0.4978s + 0.0039) 



from the Inverse Laplace Transform (Section 6.3). Using numde, sym2poly and residue, see fib4.m for 
details, returns 

/ 0.0029 \ 

262.8394 

-474.1929 

."i = -1.0857 

-9.0930 

-0.3326 

\ 211.3507 ) 



3 This content is available online at <http://cnx.Org/content/ml0523/2.3/>. 
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and 

/ -1.3565 \ 

-0.2885 

-0.1667 

l>\. = -0.1667 

-0.1667 

-0.1667 

\ -0.0100 / 

You will be asked in the exercises to show that this indeed jibes with the 



Xi{t) 



— 329£ / 

e^r I 262.842cosh 



211.35eiM - (0.0554t 3 + 4.5464t 2 + 1.085t + 474.19) e^ 
™^ l -262.836sinh('> / ^* N 



16 



16 



achieved in the Laplace Transform (Section 6.2) via ilaplace. 

8.4 Exercises: Complex Integration 4 

1. Let us confirm the representation of this Cauchy's Theorem equation (8.7) in the matrix case. More 
precisely, if $ (2) = (zl — B)~ is the transfer function associated with B then this Cauchy's Theorem 
equation (8.7) states that 

h dj 
j=lk=l \ Z ~ A 3) 



where 



<I> 



3,k 



1 

2~7ri 



$ (z) (z - X 3 ) k X dz 



.15) 



Compute the §jj, per (8.15) for the B in this equation (7.36) from the discussion of Complex Dif- 
ferentiation. Confirm that they agree with those appearing in this equation (7.40) from the Complex 
Differentiation discussion. 

2. Use this inverse Laplace Transform equation (8.14) to compute the inverse Laplace transform of 

1 

s 2 +2s+2 - 

3. Use the result of the previous exercise to solve, via the Laplace transform, the differential equation 



dt 



(x)(i) + x(i) = e~ t sin(t), x(0)=0 



Hint: Take the Laplace transform of each side. 

4. Explain how one gets from t\ and p\ to x\ (t). 

5. Compute, as in f ib4.m, the residues of C (X2 (s)) and C (x^ (s)) and confirm that they give rise to the 
X2 it) and X3 it) you derived in the discussion of Chapter 1 (Section 2.1). 



4 This content is available online at <http://cnx.Org/content/ml0524/2.4/>. 



Chapter 9 

The Eigenvalue Problem 



9.1 Introduction 1 

9.1.1 Introduction 

Harking back to our previous discussion of The Laplace Transform (Section 6.2) we labeled the complex 
number A an eigenvalue of B if XI — B was not invertible. In order to find such A one has only to find 
those s for which (si — B)~ is not defined. To take a concrete example we note that if 



B 



(l 



U 



then 



(sI-BY 



(*■ 



i)V 



/(« 



V 



1 
1 



l)(s 







2/ 



(9.1) 



s-2 
(s-l)(s- 




(*■ 



(9.2) 



I) 2 ) 



and so Ai = 1 and A2 = 2 are the two eigenvalues of B. Now, to say that Xj I — B is not invertible is to 
say that its columns are linearly dependent, or, equivalently, that the null space JV (Xjl — B) contains more 
than just the zero vector. We call JV (Xjl — B) the jth eigenspace and call each of its nonzero members a 
jth eigenvector. The dimension of JV (Xjl — B) is referred to as the geometric multiplicity of Xj. With 
respect to B above, we compute JV (X\I — B) by solving (I — B) x = 0, i.e., 



/ 


-1 




( Xl 1 




( ° \ 










x 2 


= 





V° 


1) 




Us J 




U/ 



Clearly 



B) 



jV(XiI 
Arguing along the same lines we also find 

,yV(X 2 I- B 



{< 



1 



1 



) T ' 



c£ 



c£ 



1 This content is available online at <http://cnx.Org/content/ml0405/2.4/>. 
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That B is 3x3 but possesses only 2 linearly eigenvectors leads us to speak of B as defective. The cause of 
its defect is most likely the fact that Ai is a double pole of (si — B) In order to flesh out that remark and 
uncover the missing eigenvector we must take a much closer look at the transfer function 

R(s) = {sI-By l 

In the mathematical literature this quantity is typically referred to as the Resolvent of B. 

9.2 The Resolvent 2 

9.2.1 The Transfer Function 

One means by which to come to grips with R (s) is to treat it as the matrix analog of the scalar function 

1 (9-3) 



s-b 



This function is a scaled version of the even simpler function j±-. This latter function satisfies the identity 
(just multiply across by 1 — z to check it) 

1 z n 

l + ^ + 2 + ... + ^"- 1 + - (9.4) 



\- z \- z 

for each positive integer n. Furthermore, if \z\ < 1 then z n — » as n — » oo and so (9.4) becomes, in the 
limit, 



1 oo 

1- Z ^ 



n=0 



the familiar geometric series. Returning to (9.3) we write 

1 \ 1 b b 11 - 1 b n 1 



s — b 1 — - s s 2 s n s n s — b 

s 

and hence, so long as \s\ > \b\ we find, 

1 1 °° 

£ 



s — b s *-^ V s 

n=0 x 

This same line of reasoning may be applied in the matrix case. That is, 

(--B)- 1 -.-'(/-f)" 1 (i + | + ... + ^ + ^-B)- 1 ) (-» 

and hence, so long as \s\ >\\ B \\ where || B \\ is the magnitude of the largest element of B, we find 

oo / r>\n 

(si - By 1 = s-i Y, (-) (9-6) 

n=0 ^ ' 

Although (9.6) is indeed a formula for the transfer function you may, regarding computation, not find it 
any more attractive than the Gauss- Jordan method. We view (9.6) however as an analytical rather than 
computational tool. More precisely, it facilitates the computation of integrals of R (s). For example, if C p is 
the circle of radius p centered at the origin and p >\\ B || then 



J (sI-By'ds = EZoB-fcJ-^ds 



(9.7) 
2-kiI 



2 This content is available online at <http://cnx.Org/content/ml0490/2.3/>. 
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This result is essential to our study of the eigenvalue problem. As are the two resolvent identities. Regarding 
the first we deduce from the simple observation 



(s 2 I - By 1 - (ail -B)~ l = {s 2 I - By 1 («! J -B-S2I+B) ( Sl I - B) 



that 



R{s 2 ) - iJ(ai) = (si - s 2 ) R (s 2 ) R (a{) 
The second identity is simply a rewriting of 

(si - B) {si - By 1 = {si - By 1 {si - B) = I 

namely, 

BR{s) = R{s)B 

= sR{s)-I 



(9i 



(9.9) 



9.3 The Partial Fraction Expansion of the Resolvent 3 

9.3.1 Partial Fraction Expansion of the Transfer Function 

The Gauss- Jordan method informs us that R will be a matrix of rational functions with a common denom- 
inator. In keeping with the notation of the previous chapters, we assume the denominator to have the h 
distinct roots, { Xj \ j = {1, . . . , h} } with associated multiplicities { rrij \ j = {1, . . . , h} }. 
Now, assembling the partial fraction expansions of each element of R we arrive at 



h rrij 



^) = E£ 



R 



j,k 



3 = \ k=l ( S ~ X j) 

where, recalling the equation from Cauchy's Theorem (8.7), the matrix Rj^ equals the following: 



R 



j,k 



1 

2^7 



R (z) {z — Xj) dz 



(9.10) 



(9.11) 



Example 9.1: Concrete Example 

As we look at this example with respect to the eigenvalue problem eqnl (9.1) and eqn2 (9.2), we 
find 

/ioo\ /oio\ /ooo\ 



«i. 



1 Ki,2 =000 and H 2 ,i = 000 

\ / \ / \0 1/ 

One notes immediately that these matrices enjoy some amazing properties. For example 



Ri, 



Ri, 



R 2 



R 2 



R\ \R 2 1 = 0, and R 2 







(9.12) 



Below we will now show that this is no accident. As a consequence of (9.11) and the first resolvent 
identity, we shall find that these results are true in general. 

Proposition 9.1: 



R 



JA 



Rji as seen above in (9.12). 



3 This content is available online at <http://cnx.Org/content/ml0491/2.5/>. 
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Proof: 

Recall that the Cj appearing in (9.11) is any circle about Xj that neither touches nor encircles any 
other root. Suppose that Cj and C/ are two such circles and C/ encloses Cj. Now 

Rn = R(z)dz = / R (z) dz 



and so 



1 



Rj.i = g / R{z)dz R (w) dw 

Ri i 2 = ~ / R(z)R (w) dwdz 

3 ' {2nif J J 

„ 2 _JL_ f f R(z)-R(w) J 

Rj i = t; I I dwdz 

(2tH) J J w-z 

Ri i = w / R(z) / dwdz— / R(w) / dzdw 

3 " {2TTif \J J W-Z J J W-Z 

2 ! f 

Rj i = / R (z) dz = Rj i 

2m J 

We used the first resolvent identity, This Transfer Function eqn (9.8), in moving from the second 
to the third line. In moving from the fourth to the fifth we used only 

dw = 2iri (9.13) 

w — z 

and 

dz = 



The latter integrates to zero because Cj does not encircle w. 

From the definition of orthogonal projections (Definition: "orthogonal projection", p. 51), which 
states that matrices that equal their squares are projections, we adopt the abbreviation 

p j = R j,i 

With respect to the product PjPk, for j ^ k, the calculation runs along the same lines. The 
difference comes in (9.13) where, as Cj lies completely outside of Cfc, both integrals are zero. 
Hence, 

Proposition 9.2: 

If j ^ k then PjP k = 0. 

Along the same lines we define 

Dj = Rj,2 

and prove 

Proposition 9.3: 

If 1 < k < rrij - 1 then £>/ = R Jik+1 . D 3 mj = 0. 
Proof: 
For k and I greater than or equal to one, 

Rj t k+iRj,i+i = 2 / R( z ) ( z ~ ^j) d z / R{ w ) ( w ~ Aj) dw 

(2tti) J J 
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R , , R , , — / / R (~\ R (iii\(" \ ) k (v \ \ l dird~ 




(2iri) J J 
1 f f R(z) - R(w) , skl A 

R , i R , i // ~ \ \ (in \ \dird~ 




-ftj,fc+l-ftj,i+l 2 / / „„ „ \" A 3> \ W A l> awa " 

(2iri) J J w — z 
1 f R(~)(~ \ ) k f ^ W ~ j ' dmd~ 1 f R(vA(v \ ) l 


Hz-xtf 


(2nifJ W( ~ J ' J w-z (2nifJ {H j) J 


w — z 



Rj,k+iRj,i+i = 72 / R(z)(z — Xj) j — dwdz — -^ / R (w) (w — Xj) / — dzdw 

Rj,k+\Rj.l+l = n — : / R{z) (z — Xj) dz = Rju+l+1 

2m J 
because 

(w ~ Aj) dw = 2ni{z-X 1 ) 1 (9.14) 

j w — z 

and 

^ ] -^-dz = 

w — z 

With k = I = 1 we have shown Rj t % = Rj,3, i- e -, Dj 2 = Rj^. Similarly, with k = 1 and / = 2 we 
find Rj^Rj,3 = Rj,Ai J - e -> Dj = Rja- Continuing in this fashion we find Rj t kRj,k+i = Rj,k+2 = j, 
or Dj +1 = Rj.k+2- Finally, at k = rrij — 1 this becomes 

DM = Rj, mj+1 = -L f R(z)(z- Xj) m >dz = 

Z7T1 J 

by Cauchy's Theorem. 

With this we now have the sought after expansion 

h 1 mj-l 

R (*) = £ t^ p > + D 7 — r^T^ fe (9 - 15) 

j=l Z ~ A ' fc=l (*~ A j) 

along with the verification of a number of the properties laid out in Complex Integration eqns 1-3 
(8.1). 



9.4 The Spectral Representation 4 

With just a little bit more work we shall arrive at a similar expansion for B itself. We begin by applying the 

second resolvent identity (9.9) to Pj. More precisely, we note that the second resolvent identity (9.9) implies 

that 

BP,, = PiB 

(9.16) 



= J c zR (z) — Idz 

P 3 B = / zR{z)dz 
Jc, 



PjB= / R(z)(z-Xj)dz + Xj j R{z)dz 



4 This content is available online at <http://cnx.Org/content/ml0492/2.3/>. 
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PiB = Dj + XjPj 



Summing this over j we find 



h 



B E^ = Ev^ + E^ ( 9 - 17 ) 

We can go one step further, namely the evaluation of the first sum. This stems from the eqn in the discussion 
of the transfer function (9.7) where we integrated R(s) over a circle C p where p >|| B ||. The connection to 
the Pj is made by the residue theorem. More precisely, 



L 



R{z)dz = 2™^ P. 

3 = 1 



h 

j 



Comparing this to the eqn (9.7) from the discussion of the transfer function we find 

h 

E^ = / ( 9 - 18 ) 

i=i 

and so (9.17) takes the form 

h h 

B = Y, X 3 P 3+T, D 3 ( 9 - 19 ) 

3=1 3=1 

It is this formula that we refer to as the Spectral Representation of B. To the numerous connections 
between the Pj and Dj we wish to add one more. We first write (9.16) as 

{B-\ 3 I)Pj=D 3 

and then raise each side to the nij power. As P" h = Pj and D™ 3 = we find 

(B - \jl) mj Pj = (9.20) 

For this reason we call the range of Pj the jth generalized eigenspace, call each of its nonzero members 
a jth generalized eigenvector and refer to the dimension of IZ(Pj) as the algebraic multiplicity of 
\j. With regard to the first example (p. 95) from the discussion of the eigenvalue problem, we note that 
although it has only two linearly independent eigenvectors the span of the associated generalized eigenspaces 
indeed fills out K 3 . One may view this as a consequence of Pi + P2 = I, or, perhaps more concretely, as 

appending the generalized first eigenvector 1 ) to the original two eigenvectors ( 1 ) 

and 1 ) . In still other words, the algebraic multiplicities sum to the ambient dimension (here 3), 
while the sum of geometric multiplicities falls short (here 2). 

9.5 The Eigenvalue Problem: Examples 5 

We take a look back at our previous examples in light of the results of two previous sections The Spectral 
Representation (Section 9.4) and The Partial Fraction Expansion of the Transfer Function (Section 9.3). 
With respect to the rotation matrix 

. 1 
B 

-1 



5 This content is available online at <http://cnx.Org/content/ml0493/2.4/>. 
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we recall, see Cauchy's Theorem eqn6 (8.6), that 

1 I « I 




(9.21) 



and so 



From mi = f«2 = 1 it follows that TZ (Pi) and 7?.(P 2 ) are actual (as opposed to generalized) eigenspaces. 
These column spaces are easily determined. In particular, 1Z (Pi) is the span of 




ei = 

while 1Z (P 2 ) is the span of 

e2 = 

To recapitulate, from partial fraction expansion one can read off the projections from which one can read off 
the eigenvectors. The reverse direction, producing projections from eigenvectors, is equally worthwhile. We 
laid the groundwork for this step in the discussion of Least Squares (Section 5.1). In particular, this Least 
Squares projection equation (5.14) stipulates that 

ei ei) ei and P 2 = e 2 (e 2 e 2 ) e 2 

As ei T ei = ei T ei = these formulas can not possibly be correct. Returning to the Least Squares discussion 
we realize that it was, perhaps implicitly, assumed that all quantities were real. At root is the notion of the 
length of a complex vector. It is not the square root of the sum of squares of its components but rather 
the square root of the sum of squares of the magnitudes of its components. That is, recalling that the 
magnitude of a complex quantity z is v^, 

(||ei||) 2 /ei T ei rather (|| e, ||) 2 / ef d 

Yes, we have had this discussion before, recall complex numbers, vectors, and matrices (Section 7.1). The 
upshot of all of this is that, when dealing with complex vectors and matrices, one should conjugate before 
every transpose. Matlab (of course) does this automatically, i.e., the ' symbol conjugates and transposes 
simultaneously. We use x H to denote 'conjugate transpose', i.e., 

x H = x T 

All this suggests that the desired projections are more likely 

P 1 = ei(ei ff ei)" 1 ei ff and P 2 = e 2 (e 2 ff e 2 )" 1 e 2 ff (9.22) 

Please check that (9.22) indeed jives with (9.21). 



102 



CHAPTER 9. THE EIGENVALUE PROBLEM 



9.6 The Eigenvalue Problem: Exercises 6 



9.6.1 Exercises 

1. Argue as in Proposition 1 (p. 97) in the discussion of the partial fraction expansion of the transfer 



function that if j ^ k then DjP^ 



PjD k 



0. 



2. Argue from this equation (9.17) from the discussion of the Spectral Representation that DjPj = 
PjD^D,. 

3. The two previous exercises come in very handy when computing powers of matrices. For example, 
suppose B is 4-by-4, that h = 2 and mi = mi = 2. Use the spectral representation of B together with 
the first two exercises to arrive at simple formulas for B 2 and B 3 . 

4. Compute the spectral representation of the circulant matrix 

/ 2 8 6 4 \ 

4 2 8 6 

6 4 2 8 

y 8 6 4 2 J 



B 



Carefully label all eigenvalues, eigenprojections and eigenvectors. 



6 This content is available online at <http://cnx.Org/content/ml0494/2.3/>. 



Chapter 10 

The Symmetric Eigenvalue Problem 



10.1 The Spectral Representation of a Symmetric Matrix 1 

10.1.1 Introduction 

Our goal is to show that if B is symmetric then 

• each Xj is real, 

• each Pj is symmetric and 

• each D, vanishes. 



Let us begin with an example. 

Example 10.1 

The transfer function of 



B 



/ 1 1 1 \ 

1 1 1 



is 



R(a) 



R(a) 



s{s-3) 
( 2/3 -1/3 - 



(s-2 

1 

V i 



1 

s-2 
1 



1/3 \ 

1/3 2/3 -1/3 
V -1/3 "1/3 -1/3 / 
1 



1 \ 

1 
s-2 j 

/ 1/3 1/3 1/3 \ 

1/3 1/3 1/3 
V 1/3 1/3 1/3 J 



R(s) 



Ai 



Pi 



A 2 



A 



and so indeed each of the bullets holds true. With each of the Dj falling by the wayside you may also expect 
that the respective geometric and algebraic multiplicities coincide. 



lr This content is available online at <http://cnx.Org/content/ml0382/2.4/>. 
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10.1.2 The Spectral Representation 

We have amassed anecdotal evidence in support of the claim that each Dj in the spectral representation 

h h 

* = £Vi+X>i ( 1(U ) 

j=l 3=1 

is the zero matrix when B is symmetric, i.e., when B = B T , or, more generally, when B = B H where 
B H = B Matrices for which B = B H are called Hermitian. Of course real symmetric matrices are 
Hermitian. 

Taking the conjugate transpose throughout (10.1) we find 

h h 

B " = EW +H D 3 H (10-2) 

3 = 1 3 = 1 

That is, the \j are the eigenvalues of B H with corresponding projections Pj and nilpotents Dj Hence, 
if B = B H , we find on equating terms that 



A; 


= ^3 


P3- 


-P H 

r 3 


D3- 


- D H 
^3 



and 

The former states that the eigenvalues of an Hermitian matrix are real. Our main concern however is with 
the consequences of the latter. To wit, notice that for arbitrary x, 

(|| Dj m *~ x x ||) 2 = x H (D j m =- 1 ) H D j m ^- 1 x 

(|| Dj" 1 '- 1 ! ||) 2 = x H D j m '- x D j m '- x x 
(|| D 3 m *- X x ||) 2 = x H D j m ^' 2 D j m ^x 

(||D^- 1 2 ;||) 2 = 

As Dj rrij ~ 1 x = for every x it follows (recall this previous exercise (list, item 3, p. 42)) that Dj mi ~ = 0. 
Continuing in this fashion we find Dj mj ~ 2 = and so, eventually, Dj = 0. If, in addition, B is real then as 
the eigenvalues are real and all the Dj vanish, the Pj must also be real. We have now established 

Proposition 10.1: 

If B is real and symmetric then 

h 

B = Y J ^P3 (10-3) 

3=1 

where the Xj are real and the Pj are real orthogonal projections that sum to the identity and 
whose pairwise products vanish. 
Proof: 

One indication that things are simpler when using the spectral representation is 

B 100 = ]T Xj W0 Pj (10.4) 

3=1 
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As this holds for all powers it even holds for power series. As a result, 

3=1 

It is also extremely useful in attempting to solve 

Bx = b 

for x. Replacing B by its spectral representation and b by lb or, more to the point by ^ . . Pjb we 
find 

h h 

Multiplying through by P\ gives \\P\x = P\b or Pice = ^. Multiplying through by the subsequent 
Pj's gives Pj-a; = -j-. Hence, 

x = E r lP ' x (io-5) 

We clearly run in to trouble when one of the eigenvalues vanishes. This, of course, is to be expected. 
For a zero eigenvalue indicates a nontrivial null space which signifies dependencies in the columns 
of B and hence the lack of a unique solution to Bx = b. 

Another way in which (10.5) may be viewed is to note that, when B is symmetric, this previous 
equation (9.15) takes the form 



1 



z - \f 3 
Now if is not an eigenvalue we may set z = in the above and arrive at 

B-'^yPj (10-6) 

3=1 3 

Hence, the solution to Bx = b is 

h 

A 



B-H = J2±P jb 

3 = 1 J 



as in (10.5). With (10.6) we have finally reached a point where we can begin to define an inverse 
even for matrices with dependent columns, i.e., with a zero eigenvalue. We simply exclude the 
offending term in (10.6). Supposing that Aj, = we define the pseudo-inverse of B to be 

h-l 

3 = 1 3 

Let us now see whether it is deserving of its name. More precisely, when b s 1Z (B) we would expect 
that x = B + b indeed satisfies Bx = b. Well 

h-\ h-l h-\ h-l 

BB+b = BJ2 Y P i h = E Y, BP3h = ? Y XjPjb = ? ?3h 

3=1 3 3=1 3 3 = 1 J 3=1 
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It remains to argue that the latter sum really is b. We know that 

h 
b=Y J p j b , be 31(B) 

3 = 1 

The latter informs us that b ± N (B T ). As B = B T , we have, in fact, that b _L N (B). As P h is 

nothing but orthogonal projection onto N (B) it follows that Pyja = and so B (B + b) = b, that is, 

x = B + b is a solution to Bx = b. The representation (10.4) is unarguably terse and in fact is often 

written out in terms of individual eigenvectors. Let us see how this is done. Note that if x € SR (Pi) 

then x = P\x and so, 

h 

Bx = BP lX = Y^ x 3 p j p i x = AiPiz = Aix 

i.e., x is an eigenvector of B associated with Ai. Similarly, every (nonzero) vector in SR(Pj) is an 
eigenvector of B associated with Xj . 

Next let us demonstrate that each element of 5ft (Pj) is orthogonal to each element of 5ft (P&) 
when j^k. If x e 5K (Pj) and y e 3? (P fe ) then 

z T y = (PjxfPky = x T P P k y = 

With this we note that if {xj t i,Xj t 2, ■ ■ ■ , x j,nj } constitutes a basis for Sft (Pj) then in fact the union 
of such bases, 

{ x j,P I (1 5= J ' < /*) a ^d (1 < p < rij)} 

forms a linearly independent set. Notice now that this set has 

r-i 

elements. That these dimensions indeed sum to the ambient dimension, n, follows directly from 
the fact that the underlying Pj sum to the n-by-n identity matrix. We have just proven 

Proposition 10.2: 

If B is real and symmetric and n-by-n, then B has a set of n linearly independent eigenvectors. 
Proof: 

Getting back to a more concrete version of (10.4) we now assemble matrices from the individual 
bases 

tij = \Xj t l, Xj,2, ■ ■ ■ ; x j,n j / 

and note, once again, that Pj = Ej(Ej Ej) Ej , and so 



B =Zl X 3 E 3 ( E 3 TE j) E . 



T 
3 



3 = 1 



I understand that you may feel a little overwhelmed with this formula. If we work a bit harder we can 
remove the presence of the annoying inverse. What I mean is that it is possible to choose a basis for 
each SR (Pj) for which each of the corresponding Ej satisfy Ej Ej = I As this construction is fairly 
general let us devote a separate section to it (see Gram-Schmidt Orthogonalization (Section 10.2)). 
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10.2 Gram-Schmidt Orthogonalization 2 

Suppose that M is an m-dimensional subspace with basis 

\X\ , . . . , x va y 
We transform this into an orthonormal basis 

{qi,...,q m } 
for M via 

1. Set yi = x\ and q 1 = j^ 

2. y 2 = x 2 minus the projection of x 2 onto the line spanned by q\. That is 

2/2 = x 2 - qi (qi T qi) qi T x 2 = x 2 - qiq\ T x 2 

Set q2 = ilflll and ^ 2 = ^1,92]- 

3. y 3 = x 3 minus the projection of x 3 onto the plane spanned by q\ and q 2 . That is 

2/3 = aj 3 - Q 2 (Q 2 T Q 2 ) Q 2 T x 3 = x 3 - q 1 q 1 T x 3 
Set q 3 = ip^TT and Q 3 = {qi, q 2 , q 3 }. Continue in this fashion through step (ra). 

• i m ) Vm = x m minus its projection onto the subspace spanned by the columns of Q m _i. That is 



2/m — X Tl 



(°6m—l\^m—l ^cm~l) ^m-1 



771—1 

XmXra / 

3=1 



T 

qj qj Xyn 



Set q n 



,, Vm u To take a simple example, let us orthogonalize the following basis for 



X\ 





W 



(l\ (l\ 



x 2 



1 

\oy 



x 3 



1 
V 1 J 



i- qi = yi = xi- 

2. V2 = x 2 - qiqi T x 2 = (o 1 J and so, q 2 = y 2 . 
3- 2/3 = x 3 - qiqi T x 3 = ( 1 J and so, g 3 = y 3 . 
We have arrived at 



?i 



1 




f ° 




f ° 





92 = 


1 


?3 = 





v y 




v°y 




V 1/ 



. Once the idea is grasped the actual calculations are best left to a machine. Matlab accomplishes this via 
the orth command. Its implementation is a bit more sophisticated than a blind run through our steps (1) 
through (to). As a result, there is no guarantee that it will return the same basis. For example 



2 This content is available online at <http://cnx.Org/content/ml0509/2.4/>. 
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>X=[1 1 1;0 1 1 ;0 1] ; 
>Q=orth(X) 

Q= 

0.7370 -0.5910 0.3280 

0.5910 0.3280 -0.7370 

0.3280 0.7370 0.5910 

This ambiguity does not bother us, for one orthogonal basis is as good as another. Let us put this into 
practice, via (10.8). 

10.3 The Diagonalization of a Symmetric Matrix 3 

By choosing an orthogonal basis { qjj, \ 1 < k < rij } for each R (Pj) and collecting the basis vectors in 

Qj = ( qj,i qj,2 ■ ■ ■ q 3 , nj ) 

we find that 

rij 

Pj = QjQj T = 2J Qi,k<lj,k T 
fc=l 

As a result, the spectral representation (Section 10.1) takes the form 



B — X)j=i ^jQjQj 

Eh , s~^ n 3 T 

,=1 A i l-,k=\ IjMj-k 



(10.7) 



This is the spectral representation in perhaps its most detailed dress. There exists, however, still another 
form! It is a form that you are likely to see in future engineering courses and is achieved by assembling the 
Qj into a single n-by-n orthonormal matrix 



Q={Qi ... Q h 

Having orthonormal columns it follows that Q T Q = I. Q being square, it follows in addition that Q T = Q^ 1 . 
Now 

Bqj,k = ^jqj,k 

may be encoded in matrix terms via 

BQ = QK (10.8) 

where A is the n-by- n diagonal matrix whose first n\ diagonal terms are Ai, whose next n^ diagonal terms 
are A2, and so on. That is, each A^ is repeated according to its multiplicity. Multiplying each side of (10.8), 
from the right, by Q T we arrive at 

B = QAQ T (10.9) 



3 This content is available online at <http://cnx.Org/content/ml0558/2.4/>. 
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Because one may just as easily write 

one says that Q diagonalizes B. 
Let us return the our example 



Q 1 BQ = A 



/ 1 1 1 \ 



B= 111 

V i i 1/ 

of the last chapter. Recall that the eigenspace associated with Ai = had 



ei,i 



and 



( ~l\ 

1 



ei,2 





V i / 

for a basis. Via Gram-Schmidt we may replace this with 



<?i,i 



1 

72 



1 



and 



V o 7 



91,2 



V6 



V 2 J 

Normalizing the vector associated with A2 = 3 we arrive at 



92,1 



V3 



1 
\ 1/ 



and hence 



<9 = ( Qi Qi 12 



1 
\/6 



/ -v/3 -1 \/2 \ 

V3 -1 \/2 
\ 2 V2 J 



and 



A 



/ \ 


\ 3 / 



(10.10) 
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Chapter 11 

The Matrix Exponential 



11.1 Overview 1 

The matrix exponential is a powerful means for representing the solution to n linear, constant coefficient, 
differential equations. The initial value problem for such a system may be written 

x' (t) = Ax(t) 

x(0) = x 
where A is the n-by-n matrix of coefficients. By analogy to the 1-by-l case we might expect 

x(t) = e At u (11.1) 

to hold. Our expectations are granted if we properly define e At . Do you see why simply exponentiating 
each element of At will no suffice? 

There are at least 4 distinct (but of course equivalent) approaches to properly defining e At . The first 
two are natural analogs of the single variable case while the latter two make use of heavier matrix algebra 
machinery. 

1. The Matrix Exponential as a Limit of Powers (Section 11.2) 

2. The Matrix Exponential as a Sum of Powers (Section 11.3) 

3. The Matrix Exponential via the Laplace Transform (Section 11.4) 

4. The Matrix Exponential via Eigenvalues and Eigenvectors (Section 11.5) 

Please visit each of these modules to see the definition and a number of examples. 

For a concrete application of these methods to a real dynamical system, please visit the Mass-Spring- 
Damper module (Section 11.6). 

Regardless of the approach, the matrix exponential may be shown to obey the 3 lovely properties 

1. A ( e At) = Ae At = e At A 

2. e A (* 1+ * 2 ) = e Atl e At2 

3. e At is nonsingular and (e At ) = e~( At ^ 

Let us confirm each of these on the suite of examples used in the submodules. 



1 This content is available online at <http://cnx.Org/content/ml0677/2.9/>. 



Ill 



112 CHAPTER 11. THE MATRIX EXPONENTIAL 

Example 11.1 
If 



then 




1- i (e M ) 



1 \ / e* 
2 / I e 2 * 




e 2 * 1 / \ e 2 * 2 



3. (e At ) 



Example 11.2 
If 



then 



„At 



, 1 
A. 

-1 



cos (£) sin (t) 



e 

-sin (t) cos (t) 



1. A {e M } = ( - sin (*) cos O \ and Ae M = ( -3- («) cos(t) 

\ —cos (t) —sin (£) / \ —cos (£) —sin (t) 

2. You will recognize this statement as a basic trig identity 

cos(ti + t 2 ) sin(ii+i 2 ) \ I cos (ii) sin(ii) \ / cos(t 2 ) sin(i 2 ) 
— sin(ii + t 2 ) cos(ti + t 2 ) I \ — sin (ti) cos (ti) J I — sin (t 2 ) cos(£ 2 ) 

3. Using the cofactor expansion we find 

A i = / cos(i) -sin(i) \ = l cos(-t) -sin (-i) \ = {At ) 
\ sin (i) cos (t) I \ sin (— t) cos(-t) 



Example 11.3 
If 



1 
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then 

1 t 
1 



e At 



1. £(e*)=| ° J |=A^ 

.-, , ' 'i - L l \ ( I li \ ( i. 1-2 

1 / V 1 / \ 1 

3. | l 'Y =( l "' |=e-(-) 
1/ \01 



11.2 The Matrix Exponential as a Limit of Powers 2 

You may recall from Calculus that for any numbers a and t one may achieve e at via 



fc— >oo V fc 

The natural matrix definition is therefore 



e at = limit (l+ — ) (11.2) 



fc— >oo V fc 

where / is the n-by-n identity matrix. 

Example 11.4 

The easiest case is the diagonal case, e.g., 

\ 2 

for then 

' M\ k = ( (! + £)* 

*; \ o (i+f) ; 

and so (recalling (11.2) above) 

y e 2t 
Note that this is NOT the exponential of each element of A. 

Example 11.5 

As a concrete example let us suppose 



e At = limit ( 7+4^1 ( n - 3 ) 



. 1 
A 

-1 



2 This content is available online at <http://cnx.Org/content/ml0683/2.7/>. 
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From 

, + A >=(- t l) 




We discern a pattern: the diagonal elements are equal even polynomials while the off diagonal 
elements are equal but opposite odd polynomials. The degree of the polynomial will grow with k 
and in the limit we 'recognize' 

At = | cos (*) sin (*) 

-sin (t) cos (t) 



Example 11.6 
If 



then 



for each value of k and so 



I - 



1 




At\ k lit 



k 



1 




11.3 The Matrix Exponential as a Sum of Powers 3 

You may recall from Calculus that for any numbers a and t one may achieve e at via 



fe=0 



3 This content is available online at <http://cnx.Org/content/ml0678/2.15/>. 
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The natural matrix definition is therefore 



00 (At)" 



' ■■ -E^ <"- 5 > 



^* - 

fc! 



fc=0 

where A = I is the n-by-n identity matrix. 

Example 11.7 

The easiest case is the diagonal case, e.g., 



for then 



and so (recalling (11.4) above) 



A 



(At) k 



1 
2 



t k 
{2t) k 




Note that this is NOT the exponential of each element of A. 

Example 11.8 

As a second example let us suppose 



A 



We recognize that its powers cycle, i.e.. 



A 



2 



A 3 



1 
-1 



-1 
-1 

-1 

1 

A^l 1 ° 

V ° l 

A*=l ° l 

\-l 

and so 

At _ 1 ! - 2T + Ir - • • • l - §r + It " ■ • • \ _ [ cos (*) sin (*) 



2! 


+ 


t 4 

4! 


( :i 




r' 


3! 




5! 



e 

- i+ |r-|T + --- X - It + £-••• / V " sin (*) cos (*) 
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Example 11.9 
If 



then 



and so 



. 1 
A- 





,,01 
A 2 = A 3 = A k 





, 1 t 

e At = I + tA : 

1 



11.4 The Matrix Exponential via the Laplace Transform 4 

You may recall from the Laplace Transform module 5 that may achieve e at via 



s — a 
The natural matrix definition is therefore 

e At = C 

where / is the n-by-n identity matrix. 

Example 11.10 

The easiest case is the diagonal case, e.g., 

1 
2 



,i -jC- 1 (^— ) (11.6) 



i( (s /_ A) -i) (n. 7) 



for then 



and so (recalling (11.6) above) 



(si -Ay' - ' - ' ° 

' ^2 



t- 1 (lh) e- 



Example 11.11 

As a second example let us suppose 



1 
-1 



and compute, in matlab, 



4 This content is available online at <http://cnx.Org/content/ml0679/2.10/>. 
5 "The Laplace Transform" <http://cnx.org/content/ml0731/latest/> 
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> inv(s*eye(2)-A) 

ans = [ s/(s~2+l), l/(s-2+l)] 
[-l/(s-2+l), s/(s~2+l)] 

3> ilaplace(ans) 

ans = [ cos(t), sin(t)] 
[-sin(t) , cos(t)] 



Example 11.12 
If 



then 



» inv(s*eye(2)-A) 

ans = [ 1/s, l/s~2] 
[ 0, 1/s] 

3> ilaplace(ans) 

ans = [ 1 , t] 
[ 0, 1] 



1 




11.5 The Matrix Exponential via Eigenvalues and Eigenvectors 6 

In this module we exploit the fact that the matrix exponential of a diagonal matrix is the diagonal matrix 
of element exponentials. In order to exploit it we need to recall that all matrices are almost diagonalizable. 
Let us begin with the clean case: if A is n-by-n and has n distinct eigenvalues, Xj, and therefore n linear 
eigenvectors, Sj, then we note that 

Asj = XjSj , j e {!,..., n} 



6 This content is available online at <http://cnx.Org/content/ml0680/2.13/>. 
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may be written 

A = SAS- 1 (11.8) 

where S = ( s\ s 2 ... s n ) is the full matrix of eigenvectors and A = diag (Ai, A 2 , . . . , A„) is the 
diagonal matrix of eigenvalues. One cool reason for writing A as in (11.8) is that 

A 2 = SAS^SAS- 1 = SA 2 S- 1 

and, more generally 

A k = SA k S~ 1 

If we now plug this into the definition in The Matrix Exponential as a Sum of Powers (Section 11.3), we 
find 

e At = Se M S -l 

where e At is simply 

diag( e ^,e^,...,e^) 

Let us exercise this on our standard suite of examples. 

Example 11.13 
If 

1 
2 

then 

S = I A = A 

and so e At = e At . This was too easy! 

Example 11.14 

As a second example let us suppose 

1 

\ _1 ° 
and compute, in matlab, 

> [S, Lam] = eig(A) 

S = 0.7071 0.7071 

+ 0.707H - 0.707H 



Lam = + l.OOOOi 

- l.OOOOi 



> Si = inv(S) 

Si = 0.7071 - 0.7071i 

0.7071 + 0.7071i 
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1 




3> simple (S*diag(exp(diag (Lam) *t)) *Si) 

ans = [ cos(t), sin(t)] 
[-sin(t) , cos(t)] 



Example 11.15 
If 



then matlab delivers 



> [S, Lam] = eig(A) 

S = 1.0000 -1.0000 
0.0000 

Lam =0 





So zero is a double eigenvalue with but one eigenvector. Hence S is not invertible and we can not 
invoke ((11.8)). The generalization of ((11.8)) is often called the Jordan Canonical Form or the 
Spectral Representation (Section 9.4). The latter reads 

ft 

i=i 
where the Aj are the distinct eigenvalues of A while, in terms of the resolvent R (z) = (zl — A) , 

1 



is the associated eigen-projection and 



, R(z)dz 
2ni ' 



D 3 = ^-j R ( z )( z - X j) dz 

is the associated eigen-nilpotent. In each case, Cj is a small circle enclosing only \j. 
Conversely we express the resolvent 

h . rrij-1 



n -. J 1 

,=i z A J i-i (z- A,-) 



i=i 



where 

rrij = dim (1Z (Pj)) 
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with this preparation we recall Cauchy's integral formula (Section 8.2) for a smooth function / 

1 /■/(*) 



/(«) 



-dz 



2iri J z — a 

where C (a) is a curve enclosing the point a. The natural matrix analog is 

-1 



f(A) = — / f(z)R(z)dz 

2m J 

where C (r) encloses ALL of the eigenvalues of A. For / (z) = e zt we find 



-,At 



I> V (^+E^ 



fc=i 



kV 



with regard to our example we find, h = 1, Ai = 0, P\ = I, m\ = 2, Di = A so 

e A * = / + tA 
Let us consider a slightly bigger example, if 



then 



A 



/ 1 1 \ 

1 
\ 2 / 



> R = inv(s*eye(3)-A) 

R = [ l/(s-l), l/(s-l)-2, 0] 

[ 0, l/(s-l), 0] 

[ 0, 0, l/(s-2)] 



(11.9) 



and so Ai = 1 and A2 = 2 while 



and so m\ = 2 



and 



Pi 



Di 



P?. 



( 1 \ 

1 
\ / 

/ 1 \ 


\ / 

/ \ 


\o 1) 
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and m 2 = 1 and D 2 = 0. Hence 



e At = e t (P 1 + tD 1 ) + e M P 2 



e l te l 



\ 



e* 
\ e 2t J 



11.6 The Mass-Spring-Damper System 7 




2c 



m 



Figure 11.1: Mass, spring, damper system 



If one provides an initial displacement, xq, and velocity, vo, to the mass depicted in Figure 11.1 then one 
finds that its displacement, x (t) at time t satisfies 



d 2 x (i) dx (i) 



+ kx(t) = 



(11.10) 



x(0) = x 

x' (0) = v 

where prime denotes differentiation with respect to time. It is customary to write this single second order 
equation as a pair of first order equations. More precisely, we set 

Ml (t) = X (t) 



and note that (11.10) becomes 



u 2 (t) = x' (t) 
mu 2 ' (t) = (- (fcui (t))) - 2cu 2 (t) 



(11.11) 



7 This content is available online at <http://cnx.Org/content/ml0691/2.4/>. 
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ui (i) = u 2 (t) 
Denoting u(t) = ( U \ (f) u 2 (t) ) we write (11.11) as 



u'(t)=Au(t) , A= \ fe 2c | (11.12) 



m m 



We recall from The Matrix Exponential module that 

u{t) = e At u{0) 

We shall proceed to compute the matrix exponential along the lines of The Matrix Exponential via Eigen- 
values and Eignevectors module (Section 11.5). To begin we record the resolvent 

— 1 / 2c + mz m 

R(z) = — = 

mz 2 + 2cz+k \ _fc mz 

The eigenvalues are the roots of mz 2 + 2cz + k, namely 

(-c) - d 



\i = , d = v c 2 — mk 



-c+d 



m 

We naturally consider two cases, the first being 

1. dy^O. In this case the partial fraction expansion of R(z) yields 

-111 d-c -m \ -ill c + d m \ -1 -1 

z-Xi2d\ k c + d J z-X 2 2d\ _ k d-c J z ~ x i z - \ 2 2 

and so e At = e Xlt P\ + e X2t P 2 . If we now suppose a negligible initial velocity, i.e., vo, it follows that 

x{t) = ^ (e Al * {d-c) + e A2t (c + d)) (11.13) 

If d is real, i.e., if c 2 > mk, then both Ai and A2 are negative real numbers and x (t) decays to 
without oscillation. If, on the contrary, d is imaginary, i.e., c 2 < mk, then 



x (t) = e~^ ( cos (\d\t) + ^sin (\d\t) 

\ \d\ 



(11.14) 



and so x decays to in an oscillatory fashion. When (11.13) holds the system is said to be overdamped 
while when (11.14) governs then we speak of the system as underdamped. It remains to discuss the 
case of critical damping. 

2. d = 0. In this case Ai = A2 = —y^i, and so we need only compute Pi and D\. As there is but one 
Pj and the Pj are known to sum to the identity it follows that P\ = I. Similarly, this equation (9.16) 
dictates that 

T 1 



Di = AP X - A1P1 = A - Ai/ : 

k 
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On substitution of this into this equation (11.9) we find 



e A * = e-V 



l + t 



\-t 



(11.15) 



/ 



Under the assumption, as above, that vq = 0, we deduce from (11.15) that 



x(t) 



l + t\l — | x 
m 
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Chapter 12 

Singular Value Decomposition 



12.1 The Singular Value Decomposition 1 

The singular value decomposition is another name for the spectral representation of a rectangular matrix. 
Of course if A is m-by-n and m / n then it does not make sense to speak of the eigenvalues of A. We 
may, however, rely on the previous section to give us relevant spectral representations of the two symmetric 
matrices 



A T A 
AA T 



That these two matrices together may indeed tell us 'everything' about A can be gleaned from 

JV (A T A) = J/ {A) 

jV (AA t ) = JV (A T ) 
K (A T A) = K (A T ) 

K (AA T ) = K {A) 

You have proven the first of these in a previous exercise. The proof of the second is identical. The row and 
column space results follow from the first two via orthogonality. 

On the spectral side, we shall now see that the eigenvalues of AA T and A T A are nonnegative and that 
their nonzero eigenvalues coincide. Let us first confirm this on the adjacency matrix associated with the 
unstable swing (see figure in another module (Figure 3.2: A simple swing)) 



(12.1) 



The respective products are 



AA 1 



1 ° 


1 \ 


-10 10 


\ 1/ 


( l ° ° \ 




2 




[o 1 ) 



1 This content is available online at <http://cnx.Org/content/ml0739/2.4/>. 
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10-10 

10 
A 1 A 

-10 10 

\ 1 ) 

Analysis of the first is particularly simple. Its null space is clearly just the zero vector while Ai = 2 and 
A2 = 1 are its eigenvalues. Their geometric multiplicities are m = 1 and rii = 2. In A T A we recognize 
the S matrix from the exercise in another module 2 and recall that its eigenvalues are Ai = 2, A2 = 1, and 
A3 = with multiplicities m = 1, ni = 2, and 713 = 1. Hence, at least for this A, the eigenvalues of AA T 
and A T A are nonnegative and their nonzero eigenvalues coincide. In addition, the geometric multiplicities 
of the nonzero eigenvalues sum to 3, the rank of A. 

Proposition 12.1: 

The eigenvalues of AA T and A T A are nonnegative. Their nonzero eigenvalues, including geometric 
multiplicities, coincide. The geometric multiplicities of the nonzero eigenvalues sum to the rank of 
A. 
Proof: 

If A T Ax = Xx then x T A T Ax = Xx T x, i.e., (|| Ax ||) = A(|| x ||) and so A > 0. A similar argument 
works for AA T . 

Now suppose that Xj > and that {xj^'v.Li constitutes an orthogonal basis for the eigenspace 1Z (Pj). 
Starting from 

A T Axj,k = ^jXj,k (12-2) 

we find, on multiplying through (from the left) by A that 

i.e., Xj is an eigenvalue of AA T with eigenvector Axj .£, so long as Axj t k / 0. It follows from the first 
paragraph of this proof that || Axj^ ||= \J~Xj, which, by hypothesis, is nonzero. Hence, 



AX " h 

y jtk = -j= , 1 < k < Uj (12.3) 

V Aj 

is a collection of unit eigenvectors of AA T associated with Xj. Let us now show that these vectors are 
orthonormal for fixed j. 

Vj iVj.k = ~\~ x j %A Axj t k = Xj iXj t k = 

We have now demonstrated that if Xj > is an eigenvalue of A T A of geometric multiplicity rij. Reversing the 
argument, i.e., generating eigenvectors of A T A from those of AA T we find that the geometric multiplicities 
must indeed coincide. 

Regarding the rank statement, we discern from (12.2) that if Xj > then Xj t k € 71 (A T A). The union of 
these vectors indeed constitutes a basis for 1Z (A T A), for anything orthogonal to each of these Xj t k necessarily 
lies in the eigenspace corresponding to a zero eigenvalue, i.e., in ,JV (A T A). As 1Z (A T A) = TZ (A T ) it follows 
that dim7?. (A T A) = r and hence the rij, for Xj > 0, sum to r. 

Let us now gather together some of the separate pieces of the proof. For starters, we order the eigenvalues 
of A T A from high to low, 

Ai > A 2 > ••• > X h 
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and write 
where 



A T A = XA„X T 



(12.4) 



X = {Xi, . . . , Xh] , Xj = {Xj t i, . . . , Xj, n . } 
and A„ is the n-by-n diagonal matrix with Ai in the first m slots, A2 in the next ni slots, etc. Similarly 

AA T = Yk m Y T 



(12.5) 



where 



Y = {Y u . . . ,Y h } , Y } ■ = {y jA , . . . , y jt „. } 



and A m is the m-by-m diagonal matrix with Ai in the first m slots, A2 in the next n<i slots, etc. The y^j. 
were defined in (12.3) under the assumption that Xj > 0. If Xj = let Yj denote an orthonormal basis for 
.yV (AA T ). Finally, call 



<Jj = y/Xj 

and let E denote the m-by-n matrix diagonal matrix with a\ in the first n\ slots, ai in the next n-i slots, 
etc. Notice that 

S T S = A n (12.6) 



(12.7) 



SS T = A m 
Now recognize that (12.3) may be written 

and that this is simply the column by column rendition of 

AX = YT, 

As XX T = J we may multiply through (from the right) by X T and arrive at the singular value decom- 
position of A, 

A = YT,X T (12.8) 

Let us confirm this on the A matrix in (12.1). We have 

/ 2 \ 
10 



A, 



10 
\ J 

( -1 1 



X 



V2 



\ 



A, 



V2 

10 1 
\ V2 J 

( 2 \ 

1 

\o 1) 
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Hence, 



(0 

Y = 1 

V° 

( V2 



V o 

and so A = YY,X T says that A should coincide with 



1 






1 







1/ 
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(12.9) 
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This indeed agrees with A. It also agrees (up to sign changes on the columns of X) with what one receives 
upon typing [Y, SIG, X] = scd(A) in Matlab. 

You now ask what we get for our troubles. I express the first dividend as a proposition that looks to me 
like a quantitative version of the fundamental theorem of linear algebra. 

Proposition 12.2: 

If YT,X T is the singular value decomposition of A then 

1. The rank of A, call it r, is the number of nonzero elements in E. 

2. The first r columns of X constitute an orthonormal basis for 1Z (A T ). The n — rlast columns 
of X constitute an orthonormal basis for jY (A). 

3. The first r columns of Y constitute an orthonormal basis for 1Z (A). The m — r last columns 
of Y constitute an orthonormal basis for jY (A t ). 

Let us now 'solve' Ax = b with the help of the pseudo-inverse of A. You know the 'right' thing to do, 
namely reciprocate all of the nonzero singular values. Because m is not necessarily n we must also be careful 
with dimensions. To be precise, let S + denote the n-by-m matrix whose first n\ diagonal elements are — , 
whose next n^ diagonal elements are — and so on. In the case that U] x = 0, set the final n^ diagonal elements 
of S + to zero. Now, one defines the pseudo-inverse of A to be 

A+ = XT,+Y T 

In the case of that A is that appearing in (12.1) we find 

/ V2 
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therefore, 



1 



\0 I ) 



in agreement with what appears from pinv(A). Let us now investigate the sense in which A+ is the inverse 
of A. Suppose that 6 € W" 1 and that we wish to solve Ax = b. We suspect that A + b should be a good 
candidate. Observe by (12.4) that 



because X T X = / 
by (12.6) and (12.7) 
because E T EE + = S T 
by (12.8) 



(A T A) A+b = XA n X T XT, + Y T b 

(A T A) A+b = Xk n Y>+Y T b 

(A T A) A+b = XT 1 T T,T,+Y T b 

(A T A) A+b = XY?Y T b 



(A T A) A+b = A T b 
that is A+b satisfies the least-squares problem A T Ax = A T b. 
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A A Basis for the Column Space 

Suppose A is m-by-n. If columns {cj | J = 1, ■■■,r} are the pivot columns of A re( j then columns 
{cj | j = 1, ...,r} of A constitute a basis for Ra(A). 

A Basis for the Left Null Space 

Suppose that A T is n-by-m with pivot indices { Cj \ j = {1, . . . , r} } and free indices 
{ Cj | j = {r + 1, . . . , m} }. A basis for <yK (A T ) may be constructed of m — r vectors 
{z 1 , z 2 , . . . , z m ~ r } where z k , and only z k , possesses a nonzero in its c r +k component. 

A Basis for the Null Space 

Suppose that A is m-by-n with pivot indices { Cj \ j = {1, . . . , r} } and free indices 
{cj | j = {r + 1, . . . , n} }. A basis for Ji (A) may be constructed of n — r vectors 
{z 1 , z 2 , . . . , z n ~ r } where z k , and only z k , possesses a nonzero in its c r +k component. 

A Basis for the Row Space 

Suppose A is m-by-n. The pivot rows of A Ic d constitute a basis for 1Z~\ (^4 T )- 
A subset S of a vector space V is a subspace of V when 

1. if x and y belong to S then so does x + y. 

2. if x belongs to S and t is real then tx belong to S. 
axial resistance 



l 
R 



P*N 



B Basis 

Any linearly independent spanning set of a subspace S is called a basis of S. 

C Capacitance of a Single Compartment 

C m = 27ra— c 
capacitance of cell body 

^cb = A c foC 

Column Space 

The column space of the m-by-n matrix S is simply the span of the its columns, i.e. 
Ra(S) = {Sx | x e R n }. This is a subspace (Section 4.7.2) of 5R m . The notation Ra stands for 
range in this context. 

Complex Addition 

z\ + z 2 = xi + x 2 + i {yi + 2/2) 
Complex Conjugation 

z\ = xi- iyi 
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Complex Division 



zi _ zi Z2 _ ^i^2+t(it(2+^(a:2i)i-a:ii)2) 
z 2 ~ z 2 z^ x 2 2 +y 2 2 



Complex Multiplication 

2^2 = (ari + ij/i) (x 2 + iy 2 ) = x\x 2 - 2/12/2 + i (xiy 2 + x 2 y{) 

D Dimension 

The dimension of a subspace is the number of elements in its basis. 

E Elementary Row Operations 

1. You may swap any two rows. 

2. You may add to a row a constant multiple of another row. 

F force balance 

1. Equilibrium is synonymous with the fact that the net force acting on each mass must vanish. 

2. In symbols, 

yi - yi - h = 
2/2-2/3-/2 = 



2/3-2/4-/3 = 



3. or, in matrix terms, By = f where 



/ = 



( h\ 

h 
\h ) 



and B 



/ 1 -1 \ 

1-10 
\ 1 -1 J 



H Hooke's Law 

1. The restoring force in a spring is proportional to its elongation. We call the constant of 
proportionality the stiffness, kj, of the spring, and denote the restoring force by y^. 

2. The mathematical expression of this statement is: j/j = kjej, or, 

3. in matrix terms: y = Ke where 

/ fci \ 

fc 2 

k 3 

\ k A J 



K 



I Inverse of S 



Also dubbed "S inverse" for short, the value of this matrix stems from watching what happens 
when it is applied to each side of Sx = f. Namely, 



(5a; = /) =► (S^Sx = S~ x f) =* (ix = 5" 1 /) => (x = 5~7) 



132 GLOSSARY 

Invertible, or Nonsingular Matrices 

Matrices that do have an inverse. 
Example: The matrix S that we just studied is invertible. Another simple example is 

1 

1 1 

L Left Null Space 

The Left Null Space of a matrix is the null space (Section 4.2) of its transpose, i.e., 

JT (A T ) = { y e R m | A T y = } 

Linear Independence 

A finite collection {s\, S2, ■ ■ ■ , s n } of vectors is said to be linearly independent when the only 
reals, {x\, X2, • ■ • , x n } for which x\ + X2 + ■ ■ ■ + x n = are x± = X2 = ■ ■ ■ = x n = 0. In other 
words, when the null space (Section 4.2) of the matrix whose columns are {si, s 2 , ■ ■ ■ , s n } 
contains only the zero vector. 

M Magnitude of a Complex Number 

\zi\ = ^z{z{= ^/x! 2 + y 1 2 
membrane resistance 

Rn 



2ttciJj 



N Nernst potentials 

RT [Na] RT [K] 

^ Na = ^ l0g lNa]- and Ek = ^ l0g lKl" 

where R is the gas constant, T is temperature, and F is the Faraday constant. 

Null Space 

The null space of an m-by-n matrix A is the collection of those vectors in W 1 that A maps to the 
zero vector in W 71 . More precisely, 

jV (A) = {xeR n | Ax = 0} 

O orthogonal projection 

A matrix P that satisfies P 2 = P is called a projection. A symmetric projection is called an 
orthogonal projection. 

P Pivot Column 

Each column of j4 to( j that contains a pivot is called a pivot column. 
Pivot Row 

Each nonzero row of A re( j is called a pivot row. 
Pivot 
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The first nonzero term in each row of A ro( j is called a pivot. 
poles 

Also called singularities, these are the points s at which C Xl (s) blows up. 

R Rank 

The number of pivots in a matrix is called the rank of that matrix. 
resistance of cell body 

Rcb = ^cbPm- 

Row Space 

The row space of the m-by-n matrix A is simply the span of its rows, i.e., 

Ra (A T ) = { A T y | y e R m } 

S singular matrix 

A matrix that does not have an inverse. 
Example: A simple example is: 

CO 

Span 

A finite collection {si, S2, ■ ■ ■ , s n } of vectors in the subspace S is said to span S if each element of 
S can be written as a linear combination of these vectors. That is, if for each s € S there exist n 
reals {x\, X2, ■ ■ ■ , x n } such that s = xi-si + x^s-i + • • • + x n s n . 

T The Row Reduced Form 

Given the matrix A we apply elementary row operations until each nonzero below the diagonal is 
eliminated. We refer to the resulting matrix as A ro( j. 

V Voltage-current law obeyed by a capacitor 

The current through a capacitor is proportional to the time rate of change of the potential across 
it. 
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Matrix Analysis 

Equilibria and the solution of linear and linear least squares problems. Dynamical systems and the eigenvalue 
problem with the Jordan form and Laplace transform via complex integration. 
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